Voxel-resolution myocardial blood flow analysis

ABSTRACT

A myocardial blood flow analysis scan includes incorporating a pharmacological kinetic model with the standard factor analysis model where each time activity curve is assumed to be a linear combination of factor curves. Pharmacological kinetics based factor analysis of dynamic structures (K-FADS-II) model can be applied, whereby estimating factor curves in the myocardium can be physiologically meaningful is provided. Additional optional aspects include performing a discretization to transform continuous-time K-FADS-II model into a discrete-time K-FADS-II model and application of an iterative Improved Voxel-Resolution myocardial blood flow (IV-MBF) algorithm. Where the model is applied without assumption that a right ventricle tissue curve and a left ventricle tissue curve obey a particular mathematical relationship, a least squares objective function can be applied to obtain estimates for parameters of the pharmacological kinetic model by applying a majorize-minimize optimization technique to iteratively estimate the curves.

RELATED APPLICATIONS

This application is a U.S. national phase application filed under 35U.S.C. § 371 of International Application Number PCT/US2014/059062,filed Oct. 3, 2014, designating the United States, which is acontinuation in part of U.S. application Ser. No. 14/008,021 filed Sep.27, 2013, now U.S. Pat. No. 9,149,244, which is the National Stage ofInternational Application No. PCT/US2012/031263, filed Mar. 29, 2012,which claims the benefit of U.S. Provisional application No. 61/468,765,filed Mar. 29, 2011, and this application also claims the benefit ofU.S. Provisional application No. 61/887,290, filed Oct. 4, 2013, each ofwhich is incorporated by reference in their entireties herein.

TECHNICAL FIELD

This invention relates generally to estimating myocardial blood flow(MBF) and more specifically to a positron emission tomography (PET)based method for estimating myocardial blood flow (MBF).

BACKGROUND

Imaging of blood flow through the heart and associated veins can improvediagnosis and treatment of cardiac diseases. In particular, estimationof myocardial blood flow or blood flow through heart muscular tissue canbe useful as described below.

By one approach, nuclear based medicine can be used to produce usefulmedical images. In such an approach, radioactive elements are introducedinto the bloodstream such that when the radioactive elements experiencea radioactive decay, the byproducts of that decay (often the reactionwith particles called positrons) can be sensed to produce an image ofthe area where the radioactive elements are placed. An example approachto this kind of imaging is called positron emission tomography (PET).Several radioactive elements, called positron emitting tracers, areavailable for these studies with the most common being ⁸²Rb and¹³N-Ammonia.

Currently, PET is a primary method for determining non invasive coronaryflow reserve. Coronary flow reserve can be defined as a ratio of amaximum hyperemic flow to baseline flow. In normal patients this ratiocan typically range between 3-5, which is a essentially a measure of thefunction of coronary circulation and is particularly useful in thedetection of early abnormalities due to coronary artery disease. Becausethe coronary flow reserve determination is a ratio, it is unaffected bya uniform reduction in both baseline and maximal flow.

Unfortunately, coronary flow reserve does not reflect true vasodilationcapacity. A reduction in coronary flow reserve could be caused either byincreased flow in the baseline state or by reduced maximum hyperemicflow. Factors that increase myocardial oxygen demand, for examplehypertension, increased left ventricular wall stress, increases ininotropic state, and tachycardia, can lead to an increased basal flow.Differentiating between this case and the reduced maximal hyperemic flowdue to significant coronary stenosis is difficult without absolutemyocardial blood flow measurements. Measurements of hyperemic blood flowin absolute units provide a more direct estimate of vasodilationcapacity. Accordingly, only by accurate determination of absolutemyocardial blood flow can the existence of uniform diffuse disease bedetermined.

Since the early 1990's there have been validated techniques forestimating absolute myocardial blood flow. Nevertheless, absolutemyocardial flow estimation has not been adopted for routine use in aclinic setting because of technical limitations. These limitations caninclude lack of technical expertise in a clinical setting, time taken toperform the calculations, and the lack of widely available commercialproducts to perform the calculations and display the results. On theother hand, numerous reports indicate the effect on absolute myocardialblood flow of various interventions or conditions. Yet, calculatingabsolute blood flow for clinical studies remains rare. The result isthat diagnostic decisions are usually based on relative myocardial bloodflow or relative changes in myocardial blood flow between rest andstress, often aided by a software tool that compares images to a normaldatabase.

There are at least three different kinetic models that have been used tounderstand the distribution over time of flow tracers in myocardialtissue. These works include spillover correction because of a finiteresolution of the scanner and because the myocardium is moving duringthe scan. In one known approach, factor analysis was used to obtainspillover independent time activity curves of the right ventricle (RV)and left ventricle (LV) and myocardial blood tissue. By using curvesgenerated from factor analysis, the spillover component in the model canbe eliminated in theory; however, factor analysis does not correct forthe under measurement due to the partial volume effect. Correction forthis would require the use of a contrast recovery coefficient. Methodsfor addressing the non-uniqueness problem of kinetic modeling have beenproposed. Also, kinetic modeling directly from sinograms from a dynamicsequence has been suggested.

In the following, let a_(ij) denote the activity in voxel i of frame j.In factor analysis, it is assumed that the activity is a linearcombination of K primary factor curves, where the summation coefficientsare

$\begin{matrix}{a_{ij} = {\sum\limits_{k = 1}^{K}{c_{ik}{f_{kj}.}}}} & (1)\end{matrix}$

The primary factor curves for this application are the right ventricularblood pool, the left ventricular blood pool, and the myocardial tissuecurve. The mathematical task is to find both the factors andcoefficients so that the linear combination of factor curves for everypixel in the image matches the measured curve as close as possible. Thisproblem is constrained by requiring that the tissue curves and thelinear coefficients are all positive.

Ammonia or Rubidium uptake is generally analyzed with a two or threecompartment model. The models all have a blood compartment in contactwith an extracellular free distribution compartment, which is in turn incontact with a metabolically trapped compartment. These models are mucheasier to calculate if it is assumed that the clearance from themetabolically trapped compartment is zero (or near zero) over theduration of the experiment. As a result, for accurate myocardial bloodflow modeling, several authors recommend collecting and analyzing onlytwo minutes of data. When using smooth data generated by averaging allpixels within a large range of interest, this is a reasonable approach.

While there have been significant advances in the art, further advancesare possible. For example, it is desirable to have a myocardial bloodflow analysis with greater accuracy than is presently known in the art.

SUMMARY

Generally speaking, pursuant to these various embodiments, techniquesfor estimating myocardial blood flow (MBF) in each voxel in themyocardium, and specifically to a method for estimating myocardial bloodflow in each voxel in the myocardium using a model using pharmacologicalkinetics based factor analysis of dynamic structures (K-FADS) and usinga discretization that transforms the continuous-time K-FADS model into adiscrete-time K-FADS model, then applying an iterative algorithm, suchas a Voxel-Resolution myocardial blood flow (V-MBF) algorithm.

In one approach, a myocardial blood flow analysis includes a processingdevice applying a pharmacological kinetic model to a data set stored ina storage device. The data set may be compiled from a PET scan or otherimaging approach that can monitor fluid flow in a voxel set. Forexample, the data set may be derived from an imaging technique based onmonitoring fluid based tracers in a left ventricle, a right ventricle,and myocardium of a patient or animal subject. By one aspect, thepharmacological kinetic model includes incorporating a model of changingconcentrations of bound fluid based tracers, unbound fluid basedtracers, and blood plasma fluid based tracers into a standard factoranalysis of dynamic structures model combined with a model of fluidbased tracer activity in the left ventricle as a time shifted anddispersed function of blood flow from the right ventricle. In anotherapproach, the tracer activity is modeled without assumption that a rightventricle tissue curve and a left ventricle tissue curve obey aparticular mathematical relationship. The processing device isconfigured to output a processed data set based on the application ofthe pharmacological kinetic model to the data set for providing arepresentation of blood flow in the myocardium. The processed data setmay be usable to create a visual representation, an audiorepresentation, a textural description of the myocardial blood flowusing known methods for conveying such information.

In other aspects, the processing device optionally estimates parametersof the standard factor analysis of dynamic structures model. Theestimating may be done by estimating maximum values of fluid basedtracer activity in one or both of the right ventricle or the leftventricle and modifying a corresponding signal vector value for the oneof the right ventricle or the left ventricle using the estimated maximumvalues of fluid based tracer activity. In still another approach, theestimating may be done by estimating a left ventricle time activitycurve, a right ventricle time activity curve, and a time activity curve,wherein the left ventricle time activity curve is assumed to beapproximately equal to a response of an mth-order, all-pole filterapplied to the right ventricle time activity curve and determining a setof parameters that produce a smallest least-squares error for thepharmacological kinetic model. This estimation may include for a giveninitial estimate right ventricle time activity curve and a given initialestimated left ventricle time activity curve, determining initialestimates for parameters of the pharmacological kinetic model.

Where the model is applied without assumption that a right ventricletissue curve and a left ventricle tissue curve obey a particularmathematical relationship, a least squares objective function can beapplied to obtain estimates for parameters of the pharmacologicalkinetic model. In one such approach, the least squares objectivefunction is minimized by applying a majorize-minimize optimizationtechnique to iteratively estimate the right ventricle tissue curve andthe left ventricle tissue curve. At initialization, initial estimatesfor an initial estimated right ventricle time activity curve and aninitial estimated left ventricle time activity curve are estimated anduse the initial estimates for the initial estimate right ventricle timeactivity curve and the initial estimated left ventricle time activitycurve to determine initial parameters of the pharmacological kineticmodel. Where the data is noisy, the processing device smoothingestimates for blood flow for a given voxel by applying a limiting factoror penalty factor on data from voxels located within a given distancefrom the given voxel. Optionally, the pharmacological kinetic modelincludes a semi-parametric model configured to drive time activitycurves for the right ventricle imaging activity from the fluid basedtracers and the left ventricle imaging activity from the fluid basedtracers to zero over time.

So configured, a more accurate derivation of myocardial blood flow ispossible through the application of these modeling techniques. Theadvantages include increased resolutions and incorporation of kinetics.Further, unlike most known iterative algorithms for MBF estimation, suchapproaches explicitly describe a general procedure for initializing suchalgorithms. Other features will become more apparent to persons havingordinary skill in the art from the following description and claims.

BRIEF DESCRIPTION OF THE DRAWINGS

The foregoing features, as well as other features, will become apparentwith reference to the description and figures below, in which likenumerals represent elements, and in which:

FIG. 1 illustrates a generalization of a 3-compartment model used tomodel Rubidium or Ammonia kinetics in the heart.

FIG. 2 illustrates results from a V-MBF algorithm applied to syntheticdata: (a) true, initial, and estimated RV factor, (b) comparison of trueTAC (computed analytically) and estimated TAC generated using thediscrete-time K-FADS model and estimated parameters, (c) same as (b)except for a different TAC, and (4) least-squares objective function asa function of iteration number.

FIG. 3 illustrates example results for plane 9 using the IV-MBFalgorithm where FIG. 3A illustrates MBF estimates displayed as an image,FIG. 3B illustrates a mask used to define myocardium voxels, FIG. 3Cillustrates the image obtained by averaging the PET images for the lastsix sub-scans, and FIG. 3D illustrates a histogram of the MBF estimates.

FIG. 4 illustrates example results for plane 16 using the IV-MBFalgorithm, where FIG. 4A illustrates MBF estimates displayed as animage, FIG. 4B illustrates a mask used to define myocardium voxels. FIG.4C illustrates the image obtained by averaging the PET images for thelast six sub-scans, and FIG. 4D illustrates a histogram of the MBFestimates.

FIG. 5 illustrates example results for a specific myocardium voxel,where FIG. 5A illustrates a plot of the measured interpolated TAC(solid) and corresponding estimated TAC using the IV-MBF algorithm(dashed-dotted), FIG. 5B illustrates signal components that combine toform estimated TAC: sampled RV and LV tissue curves, and curves due toactivity in the free and bound states are represented by the dashed,dotted, solid, and dashed-dotted lines, respectively, and 5C illustratesthe LS objective function.

FIG. 6 illustrates an embodiment of a illustrative system that canincorporate the present embodiments.

Skilled artisans will appreciate that elements in the figures areillustrated for simplicity and clarity and have not necessarily beendrawn to scale. For example, the dimensions and/or relative positioningof some of the elements in the figures may be exaggerated relative toother elements to help to improve understanding of various embodimentsof the present invention. Also, common but well-understood elements thatare useful or necessary in a commercially feasible embodiment are oftennot depicted in order to facilitate a less obstructed view of thesevarious embodiments. It will further be appreciated that certain actionsand/or steps may be described or depicted in a particular order ofoccurrence while those skilled in the art will understand that suchspecificity with respect to sequence is not actually required. It willalso be understood that the terms and expressions used herein have theordinary technical meaning as is accorded to such terms and expressionsby persons skilled in the technical field as set forth above exceptwhere different specific meanings have otherwise been set forth herein.

DETAILED DESCRIPTION

A method and apparatus for estimating myocardial blood flow (MBF) ineach voxel in the myocardium is described. The algorithm is based on afactor analysis of dynamic structures (FADS) model that has beenenhanced to constrain the factor analysis curves to be physiologicallyappropriate.

Turning now the figures, FIG. 1 illustrates a generalization of a3-compartment model used to model Rubidium or Ammonia kinetics in theheart. Tracer in the second compartment is freely diffusible, while thetracer in the third compartment is trapped. The unknown parameters k₁,k₂, and k₃ are first order rate constants that describe the tracermovement between the compartments. The parameter k₁ (ml/min/g) isidentified as myocardial blood flow. A first approach to analyzing MBFwill now be described.

I. A First Approach

The kinetic behavior of ammonia in the myocardium is modeled with thecompartment model shown in FIG. 1 and the following differentialequations (note: assume k₄=0):

$\begin{matrix}{\frac{d\;{C_{u}(t)}}{d\; t} = {{k_{1}{C_{p}(t)}} - {\left( {k_{2} + k_{3}} \right){C_{u}(t)}}}} & (2) \\{{\frac{d\; C_{b}}{d\; t} = {k_{3}{C_{u}(t)}}},} & (3)\end{matrix}$where C_(p) is the ammonia concentration in blood plasma, C_(u) is theconcentration of the free (i.e., unbound) ammonia, and C_(b) is theconcentration of the trapped (i.e., bound) ammonia.

In one approach, a model is applied that incorporates a pharmacologicalkinetic model with the standard FADS model, which is a model where eachtime activity curve is assumed to be a linear combination of factorcurves. The resulting model, which can be called the pharmacologicalkinetics based factor analysis of dynamic structures (K-FADS) model,provides a means for estimating factor curves in the myocardium that arephysiologically meaningful. Further, a discretization is performed totransform continuous-time K-FADS model into a discrete-time K-FADSmodel. It should be noted that there is a simple relationship betweenthe discrete-time and continuous-time K-FADS parameters.

Next, an iterative algorithm can be applied, such as theVoxel-Resolution myocardial blood flow (V-MBF) algorithm. This algorithmcan iteratively estimate the MBF for each voxel in the myocardium. TheV-MBF algorithm is reliably initialized using an input-output systemidentification method described by Steiglitz-McBride. This method can beapplicable to a class of discrete-time systems that includes thediscrete K-FADS model. This V-MBF algorithm was evaluated subjectivelyand objectively through experiments conducted using synthetic data andfound to be reliable in considered test scenarios. Accordingly, themethod disclosed herein is feasible for determining physiologicallymeaningful estimates of absolute MBF.

The initial model is a continuous-time K-FADS model from which adiscrete-time K-FADS model can be obtained by applying a bilineartransform to the continuous-time K-FADS model. It should be noted thatthere is a simple relationship between the discrete-time andcontinuous-time K-FADS parameters. Using a systems theory framework, theproblem of estimating the discrete-time K-FADS parameters is related tothe problem of identifying a discrete-time system from given input andoutput data (i.e., input-output system identification). The V-MBFalgorithm iteratively estimates the MBF for each voxel in themyocardium. The V-MBF algorithm can be initialized using an input-outputsystem identification method, such as one described bySteiglitz-McBride, which is applicable to a class of discrete-timesystems that includes the discrete K-FADS model. Each of these aspectswill be described in turn with respect to a first approach to estimatingthe MBF for voxels of the cardiac system.

Continuous-Time K-FADS Model

For convenience, the RV factor, LV factor (i.e., Cp). Cu, and Cb can bedenoted by the continuous-time functions f₁(t), f₂(t), f_(i,3)(t), andf_(i,4)(t), respectively, where i is the voxel index. With this notationit follows that f₁ and f₂ represent the activity concentration of theammonia in the right ventricle and left ventricle, respectively.Additionally, f_(i,3)(t), and f_(i,4)(t) represent the activityconcentration of free and trapped ammonia in myocardial tissue,respectively. It can be seen that it is assumed that the RV (rightventrical) and LV (left ventrical) factors are spatially constant (i.e.,RV and LV factor are voxel independent). In the continuous-time K-FADSmodel the LV factor is modeled as a time shifted and dispersed functionof the RV factorf ₂(t)=γf ₁(t)*exp(−β(t−τ))u(t−τ)  (4)where γ, β, and τ (tau) are the unknown gain, time constant, and delayof the LV, respectively. Specifically, τ (tau) accounts for the factthat the ammonia activity first appears in the right ventricle and then,after a period of time, appears in the LV. The function u is the unitstep function, and the notation * denotes the continuous-timeconvolution operator. The model for the LV factor in (4) can bemotivated by observations of “isolated” RV and LV factors obtained fromdynamic PET sequences and, in part, by the need for mathematicaltractability. For example, consider if parameters k₁, k₂, and k₃ arepixel dependent and voxel i lies in the myocardium, then, applying theLaplace transform to (3) can lead to the following expressions for theactivity concentration of the free and trapped ammonia in voxel if _(i,3)(t)=k _(i,1) f ₂(t)*exp(−(k _(i,2) +k _(i,3))t)u(t)  (5)f _(i,4)(t)=k _(i,3) f _(i,3)(t)*u(t).  (6)It is noted that k_(1,1), k_(2,1), . . . , k_(1,1) are the preferred MBFparameters. In keeping with the assumptions behind (1), the activity forthe ith pixel can be expressed asa _(i)(t)=c _(i,1) f ₁(t)+c _(i,2) f ₂(t)+c _(i,3) f _(i,3)(t)+c _(i,4)f _(i,4)(t).  (7)

The first term in (7) can be identified as the amount of spillover fromthe right ventricle, and the second term can be a combination of theammonia activity in blood vessels within the myocardium and spilloverfrom the left ventricle. More specifically, the constant c_(i,1)accounts for the amount of the measured radioactivity in voxel in thecase of a PET scan that is due to the blood plasma in the RV. Further,the constant c_(i,2) accounts for the amount of the measuredradioactivity in voxel i that is due to blood plasma in the LV (i.e., LVspill over) and blood plasma in the blood vessels of the myocardium. Forthis approach, it is assumed that c_(i,2)=0.05. The third and fourthterms in (7) are the activity of the free and trapped ammonia in themyocardial tissue, respectively. The coefficients c_(i,3) and c_(i,4)represent the fractional volume of voxel i that can be occupied by theradiotracer in either the free or trapped states. Given the free spacefor water in myocardial tissue is approximately 80 percent, it isassumed that c_(i,3)=c_(i,4)=0.8.

Within the descriptions herein, it can be convenient to let f(t)=f₁(t)and c_(i)=c_(i,1). From equations (4), (5), (6), and straightforwardcalculations, the functions f_(i,3)(t) and f_(i,4)(t) can be written as:

$\begin{matrix}{\mspace{79mu}{k_{i} = {k_{i,2} + k_{i,3}}}} & (8) \\{\mspace{79mu}{{f_{i,3}(t)} = {{- {f(t)}}*\left\lbrack {{\exp\left( {- {\beta\left( {t - \tau} \right)}} \right)} - {\exp\left( {- {k_{i}\left( {t - \tau} \right)}} \right)}} \right\rbrack{u\left( {t - \tau} \right)}}}} & (9) \\{{f_{i,4}(t)} = {\gamma\; k_{i,1}k_{i,3}{f(t)}*\left\lbrack {\frac{\exp\left( {- {\beta\left( {t - \tau} \right)}} \right)}{\beta\left( {\beta - k_{i}} \right)} - \frac{\exp\left( {- {k_{i}\left( {t - \tau} \right)}} \right)}{k_{i}\left( {\beta - k_{i}} \right)} + \frac{1}{\beta\; k_{i}}} \right\rbrack{u\left( {t - \tau} \right)}}} & (10)\end{matrix}$

By substituting these equations into (7) it can be shown that

$\begin{matrix}{{{a_{i}(t)} = {{c_{i}{f(t)}} + {{f(t)}*\left\lbrack {b_{i,1} + {b_{i,2}{\exp\left( {- {\beta\left( {t - \tau} \right)}} \right)}} + {b_{i,3}{\exp\left( {- {k_{i}\left( {t - \tau} \right)}} \right)}}} \right\rbrack{u\left( {t - \tau} \right)}}}},} & (11) \\{\mspace{79mu}{where}} & \; \\{\mspace{79mu}{b_{i,1}\overset{\Delta}{=}\frac{\gamma\; c_{i,4}k_{i,1}k_{i,3}}{\beta\; k_{i}}}} & (12) \\{\mspace{79mu}{b_{i,2}\overset{\Delta}{=}{\gamma\left( {c_{i,2} - \frac{c_{i,3}k_{i,1}}{\beta - k_{i}} + \frac{c_{i,4}k_{i,1}k_{i,3}}{\beta\left( {\beta - k_{i}} \right)}} \right)}}} & (13) \\{\mspace{79mu}{b_{i,3}\overset{\Delta}{=}{\frac{\gamma\; k_{i,1}}{\beta - k_{i}}\left( {c_{i,3} - \frac{c_{i,4}k_{i,3}}{k_{i}}} \right)}}} & (14)\end{matrix}$

It can also be shown from straightforward calculations that

$\begin{matrix}{k_{i,1} = {\frac{c_{i,2}}{c_{i,3}}{\frac{{\left( {\beta - k_{i}} \right)b_{i,3}} + {\beta\; b_{i,1}}}{b_{i,1} + b_{i,2} + b_{i,3}}.}}} & (15)\end{matrix}$Thus, given estimates for the parameters β, {k_(i)}, {c_(i)}, {b_(i,1)},{b_(i,2)}, and {b_(i,3)}, and using the assumed values for {c_(i,2)} and{c_(i,3)}, the MBF parameters {k_(i,1)} can be estimated from equations(15). As stated above, it is assumed that c_(i,2)=0.05 and c_(i,3)=0.8for all i.

Discrete-Time K-FADS Model

The next aspect of this approach is to address the problem of estimatingthe parameters β, {k_(i)}, {c_(i)}, {b_(i,1)}, {b_(i,2)}, and {b_(i,3)}from discrete-time Time Activity Curve (TAC) data becausecontinuous-time TAC data is not available in practice.

In practice, only samples of the TACs are available so let y_(i)[n],i=1, 2, . . . I and x[n] denote the discrete-time signals obtained bysampling the ith TAC a_(i)(t) and RV factor respectivelyy _(i) [n]

a _(i)(nT), n=0, 1, . . . , N−1  (16)x[n]

f(nT), n=0, 1, . . . , N−1,  (17)

where f_(s)

1/T is the sampling rate and T is the sampling interval. It is notedthat in applications where scan durations for dynamic sequence PETprotocols are not uniform, the assumption that the TACs are sampleduniformly is inappropriate. However, uniform samples of the TACs can beobtained from non-uniform samples of the TACs via a suitableinterpolation. It follows that, in order to estimate the MBF parameters,the parameters β, {k_(i)}, {c_(i)}, {b_(i,1)}, {b_(i,2)}, and {b_(i,3)}are estimated from the data y_(i)[n], n=0, 1, . . . , N−1. It can beobserved that the parameters {c_(i)} are really nuisance parametersbecause they do not show up in the expression for ki,1 (see equation(15)).

It is of interest to determine a discrete-time system with the propertythat its response to the discrete-time RV factor x[n] closelyapproximates the ith discrete-time TAC y_(i)[n]. The bilineartransformation is a way to transform a linear time-invariantcontinuous-time system into a linear time-invariant discrete-timesystem. A limitation of the bilinear transform is that a delay in acontinuous-time system must be an integer multiple of the samplinginterval. Taking the Laplace transform of (11), we get the followingrelationship

$\begin{matrix}{{{{A_{i}(s)} = {{c_{i}{F(s)}} + {{F(s)}\left( {{b_{i,1}\frac{1}{s}} + {b_{i,2}\frac{1}{s + \beta}} + {b_{i,3}\frac{1}{s + {ki}}}} \right){\exp\left( {{- s}\;\tau} \right)}}}},\mspace{20mu}{i = 1},2,\ldots\mspace{14mu},I}\;} & (18)\end{matrix}$

As a result, it follows that the system function of the overallcontinuous-time system is given by

$\begin{matrix}{{H_{i,{tot}}(s)}\overset{\bigtriangleup}{=}{\frac{A_{i}(s)}{F(s)} = {c_{i} + {\left( {{b_{i,1}\frac{1}{s}} + {b_{i,2}\frac{1}{s + \beta}} + {b_{i,3}\frac{1}{s + {ki}}}} \right){\exp\left( {{- s}\;\tau} \right)}}}}} & (19)\end{matrix}$

Assuming the delay t is a multiple of the sampling interval T, thesystem function of the desired discrete-time system, H_(i)(z), can beobtained by applying the bilinear transformation to the overallcontinuous-time system H_(i,tot)(s)

$\begin{matrix}{{H_{i}(z)} = {{c_{i} + {\left( \left. \left\lbrack {{b_{i,1}\frac{1}{s}} + {b_{i,2}\frac{1}{s + \beta}} + {b_{i,3}\frac{1}{s + {ki}}}} \right\rbrack \right|_{s = {\frac{2}{T}\frac{1 - s^{- 1}}{1 + s^{- 1}}}} \right)z^{- d}}} =}} & (20) \\{c_{i} + {\left( {{b_{i,1}^{\prime}\frac{1 + z^{- 1}}{1 - {r_{1}z^{- 1}}}} + {b_{i,2}^{\prime}\frac{1 + z^{- 1}}{1 - {r_{2}z^{- 1}}}} + {b_{i,3}^{\prime}\frac{1 + z^{- 1}}{1 - {r_{i,3}z^{- 1}}}} +} \right)z^{- d}}} & (21)\end{matrix}$where τ=dT for some integer d, r₁=1, r₂=(2/T+β)⁻¹(2/T−β),r_(i,3)=(2/T+k_(i))⁻¹(2/T−k_(i)), b′_(i,1)=b_(i,1)(2/T)⁻¹,b′_(i,2)=b_(i,2)(2/T+β)⁻¹, and b′_(i,3)=b_(i,3)(2/T+k_(i))⁻¹.

It is noted that for this approach, the z-transform is used. Here letg[n] be an arbitrary discrete-time sequence. The z-transform g[n] can bedefined as

${G(z)} = {\sum\limits_{n = {- \infty}}^{\infty}\;{{g\lbrack n\rbrack}{Z^{- n}.}}}$It should be noted that the delay term (i.e., z^(−d) term) in equation(21) follows from the τ second delay in the continuous-time K-FADSmodel, which is the delay between activity appearing in the right andleft ventricles.

As known in the art, the bilinear transformation has the property thatit maps a stable continuous-time system into a stable discrete-timesystem. Moreover, the bilinear transformation avoids the problem ofaliasing by mapping the jΩ axis into the unit circle of the complexplane. However, “frequency warping” can occur as a result of mapping theentire jΩ axis into the unit circle. Note, the frequency warping problemcan be ameliorated by choosing a sufficiently high sampling rate f_(s).

It follows from the definition in equation (21) that the discrete-timeK-FADS model for the ith TAC can be represented by the followinginput-output relationship

$\begin{matrix}{{Y_{i}(z)} = {{{X(z)}{H_{i}\left( {{z;d},r_{2},r_{i,3},\theta_{i}} \right)}} = {{c_{i}{X(z)}} + {{{X(z)}\left\lbrack {{b_{i,1}^{\prime}\frac{1 + z^{- 1}}{1 - {r_{1}z^{- 1}}}} + {b_{i,2}^{\prime}\frac{1 + z^{- 1}}{1 - {r_{2}z^{- 1}}}} + {b_{i,3}^{\prime}\frac{1 + z^{- 1}}{1 - {r_{3}z^{- 1}}}}} \right\rbrack}z^{- d}}}}} & (22)\end{matrix}$where θ_(i)

[b′_(i,1), b′_(i,2), b′_(i,3), c_(i)](recall r_(i)=1). The notationH_(i)(z; d, r₂, r_(i,3), θ_(i)) explicitly illustrates the dependence ofthe ith system function on the unknown parameters. For the discussionbelow, it is beneficial to define the following notation:θ_(r) ₃

[r_(1,3), r_(2,3), . . . , r_(1,3)]θ_(b′)

[b′_(1,1), b′_(2,1), . . . , b′_(I,1) , b′ _(1,2), b′_(2,2), . . .,b′_(I,2), b′_(1,3), b′_(2,3), . . . , b′_(I,3)]c

[c₁, c₂, . . . , c_(I)].  (23)

The problem of interest is to estimate parameters of the discrete-timeK-FADS model from the sampled TACs y_(i)[n], I=1, 2, . . . , I.Therefore, the V-MBF algorithm of this approach solves the followingleast-squares problem

$\begin{matrix}{{{(P)\left( {\hat{x},\hat{d},{\hat{r}}_{2},{\hat{\theta}}_{r_{3}},{\hat{\theta}}_{b^{\prime}},\hat{c}} \right)} = {\arg\;{\min\limits_{{x \geq 0},{\theta \in s_{\theta}}}{\phi\left( {x,d,r_{2},\theta_{r_{3}},\theta_{b^{\prime}},c} \right)}}}},} & (24) \\{where} & \; \\{{\phi\left( {x,d,r_{2},\theta_{r_{3}},\theta_{b^{\prime}},c} \right)}\overset{\bigtriangleup}{=}{\sum\limits_{i = 1}^{I}\;{\sum\limits_{n = 1}^{N}\;\left( {{y_{i}\lbrack n\rbrack} - {{x\lbrack n\rbrack}*{h_{i}\left\lbrack {{n;d},r_{2},r_{i,3},\theta_{i}} \right\rbrack}}} \right)^{2}}}} & (25) \\{\theta\overset{\bigtriangleup}{=}{\left\lbrack {d,r_{2},\theta_{r_{3}},\theta_{b^{\prime}},c} \right\rbrack.}} & (26)\end{matrix}$

Additionally, S_(θ) is a feasible set of the parameters θ and h_(i)[n;d, r₂, r_(i,3), θ_(i)] is the inverse Z-transform of H_(i)(z; d, r₂,r_(i,3), θ_(i)).

In the development of the V-MBF algorithm, the problem (P) is simplifiedby assuming that the discrete-time d delay is known. Next, thisassumption is removed by estimating the parameter d along with the otherparameters of the discrete-time K-FADS model. An initialization methodis applied that exploits the well known Steiglitz-McBride algorithm.

A. Discrete-Time Delay Known

To minimize the objective function

in (24), an algorithm is developed based on the group coordinate-descentmethod. By use of the term group coordinate-descent method it isunderstood that, in a cyclic fashion, the objective function

is minimized with respect to a set of parameters while the otherparameters are fixed.

Let do denote the known discrete-time delay. Given initial estimates,x⁽⁰⁾, r₂ ⁽⁰⁾, θ_(r) ₃ ⁽⁰⁾, θ_(b′) ⁽⁰⁾, c⁽⁰⁾ an algorithm for solving (P)proceeds as follows:

For m = 0, 1, . . . , M − 1 or repeat until some chosen criterion ismet: Step 1$x^{({m + 1})} = {\arg\mspace{11mu}{\min\limits_{x \geq 0}\;{\phi\left( {x,d_{0},r_{2}^{(m)},\theta_{r_{3}}^{(m)},\theta_{b^{\prime}}^{(m)},c^{(m)}} \right)}}}$Step 2$\left( {r_{2}^{({m + 1})},\theta_{r_{3}}^{({m + 1})}} \right) = {\arg{\min\limits_{{0 < r_{2}},{\theta_{r_{3}} < 1}}{\phi\left( {x^{({m + 1})},d_{0},r_{2},\theta_{r_{3}},\theta_{b^{\prime}}^{(m)},c^{(m)}} \right)}}}$(27) Step 3$\left( {\theta_{b^{\prime}}^{({m + 1})},c^{({m + 1})}} \right) = {\arg\mspace{11mu}{\min\limits_{0 \leq c \leq 1}\;{\phi\left( {x^{({m + 1})},d_{0},r_{2}^{({m + 1})},\theta_{r_{3}}^{({m + 1})},\theta_{b^{\prime}},c} \right)}}}$(28) end

Solution to Step 1 of the V-MBF Algorithm with Known Delay

In the solution to Step 1 of the V-MBF algorithm with known delay d, itis convenient to denote the next estimate for the RV factor as x

x^((m+1)). Given its simplicity, the coordinate descent algorithm isused to iteratively determine x. Let e₁, e₂, . . . , e_(N) be the searchdirections, where e_(j) is an N×1 vector of zeros except for a 1 at thejth position. The steps of the coordinate descent algorithm fordetermining the update x^((m+1)) are as follows: Let z₁=x^((m)), x⁽⁰⁾=x^((m)), j=1, and 1=0, and choose ∈>0

$\left. {{{{Step}\mspace{14mu} 1.1\mspace{14mu}{Determine}\mspace{14mu}\lambda_{j}} = {\arg\;{\min\limits_{\lambda}{\phi\left( {z_{j} + {\lambda\; e_{j}}} \right)}}}},d_{0},r_{2}^{(m)},\theta_{r_{3}}^{(m)},\theta_{b^{\prime}}^{(m)},c^{(m)}} \right)$

Step 1.2 Let z_(j+1)=z_(j)+λ_(j)e_(j). If jth component of z_(j+1) isnegative, then this value is set to zero. Note, this operation accountsfor the nonnegativity constraint of the discrete-time RV factor. If j<N,then increment j and repeat Step 1.1. Otherwise, if j=N, then go to Step1.3.

Step 1.3 Let x ^((l+1))=z_(N+1). If ∥x ^((l+1))−x ^((l))∥<∈ (or aspecified number of desired iterations is reached), then stop and letx^((m+1))=x ^((l+1)). Else, let z₁

x ^((l+1)) and j=1, increment l, and go to Step 1.1. Note, the step sizeλ_(j) has a closed form expression.

Solution to Step 2 of the V-MBF Algorithm with Known Delay

Again, the simplicity of the coordinate descent method is exploited tocompute a solution to the problem in Step 2. However, the coordinatedescent method is expressed in a manner that is more convenient for thisproblem:

repeat until ∥θ_(r) ^((m+1)) − θ_(r) ^((m))∥ < ϵ where θ_(r) ^((m)) 

 [r₂ ^((m)), r_(1,3) ^((m)), r_(2,3) ^((m)), . . . , r_(I,3) ^((m))]Step 2.1:$r_{2}^{({m + 1})} = {\arg\mspace{11mu}{\min\limits_{0 < r_{2} < 1}\;{\phi\left( {x^{({m + 1})},d_{0},r_{2},\theta_{r_{3}}^{(m)},\theta_{b^{\prime}}^{(m)},c^{(m)}} \right)}}}$Step 2.2: for i = 1, 2, . . . , I, (29)$r_{i,3}^{({m + 1})} = {\arg{\min\limits_{0 < r_{3} < 1}{\sum\limits_{n = 1}^{N}\left( {{y_{i}\lbrack n\rbrack} - {{x^{({m + 1})}\lbrack n\rbrack}*{h_{i}\left\lbrack {{n;d},r_{2}^{({m + 1})},r_{3},\theta_{i}^{({m + 1})}} \right\rbrack}}} \right)^{2}}}}$end end

It should be observed that in Step 2.2 advantage is taken of the factthat the objective function

is de-coupled in terms of the parameters r_(1,3), r_(2,3), . . . ,r_(I,3). Also, a 1-D line search algorithm such as the golden sectionmethod can be used to solve the 1D minimization problems in Steps 2.1and 2.2.

Solution to Step 3 of the V-MBF Algorithm with Known Delay

Referring to equation (25), it follows that the problem in equation (28)is equivalent to the following problems, i=1,2, . . . , I,

$\begin{matrix}{\theta_{i}^{({m + 1})} = {\arg_{0 \leq c_{i} \leq 1}{\min\limits_{b_{i,1}^{\prime},b_{i,2}^{\prime},b_{i,3}^{\prime}}{\sum\limits_{n = 1}^{N}\;\left( {{y_{i}\lbrack n\rbrack} - {c_{i}{x^{({m + 1})}\lbrack n\rbrack}} - {b_{i,1}^{\prime}{w_{1}^{({m + 1})}\lbrack n\rbrack}} - {b_{i,2}^{\prime}{w_{2}^{({m + 1})}\lbrack n\rbrack}} - {b_{i,3}^{\prime}{w_{3}^{({m + 1})}\lbrack n\rbrack}}} \right)^{2}}}}} & (30)\end{matrix}$where, w_(p) ^((m+1))[n], p=1, 2, is the inverse Z-transform of

$\begin{matrix}{{{W_{p}^{({m + 1})}(z)} = {{X^{({m + 1})}(z)}\frac{1 + {z^{- 1}z^{- d}}}{1 - {r_{p}z^{- 1}}}}},{p = 1},2,} & (31)\end{matrix}$and w_(i,3) ^((m+1))[n] is the inverse Z-transform of

$\begin{matrix}{{W_{i,3}^{({m + 1})}(z)} = {{X^{({m + 1})}(z)}{\frac{1 + {z^{- 1}z^{- d}}}{1 - {r_{i,3}^{({m + 1})}z^{- 1}}}.}}} & (32)\end{matrix}$

An optimization problem in equation (30) is a linear least-squaresproblem under the constraint 0≤c_(i)≤1. If the constraint on c_(i)ignored, then the update θ_(i) ^((m+1)) could be computed by solving thenormal equations associated with the least-squares objective function in(30). Thus, if the solution to an unconstrained version of theleast-squares problem in (30) is such that 0≤c_(i) ^((m+1))≤1 for all i,then no other steps should be necessary. Alternatively, if without lossof generality, the constraint is not satisfied for i=i₀, then additionalsteps should be taken. A straightforward strategy could be to firstcompute the updates b′_(i) ₀ _(,1) ^((m+1)), b′_(i) ₀ _(,2) ^((m+1)),and b′_(i) ₀ _(,3) ^((m+1)) by solving the normal equations associatedwith the following least-squares problem

$\begin{matrix}{\left( {b_{i_{0},1}^{\prime{({m + 1})}},b_{i_{0},2}^{\prime{({m - 1 - \pi})}},b_{i_{0},3}^{\prime{({m + 1})}}} \right) = {\arg\;{\min\limits_{b_{i_{0},1}^{\prime},b_{i_{0},2}^{\prime},b_{i_{0},3}^{\prime}}{\sum\limits_{n = 1}^{N}\left( {{y_{i_{0}}\lbrack n\rbrack} - {c_{i_{0}}^{(m)}{x^{({m + 1})}\lbrack n\rbrack}} - {b_{i_{0},1}^{\prime}{w_{1}^{({m + 1})}\lbrack n\rbrack}} - {b_{i_{0},2}^{\prime}{w_{2}^{({m + 1})}\lbrack n\rbrack}} - {b_{i_{0},3}^{\prime}{w_{3}^{({m + 1})}\lbrack n\rbrack}}} \right)^{2}}}}} & (33)\end{matrix}$

The update c_(i) ₀ ^((m+1)) could then be computed using the coordinatedescent method

$\begin{matrix}{c_{i_{0}}^{({m + 1})} = {\arg\;{\min\limits_{0 \leq c \leq 1}{\sum\limits_{n = 1}^{N}\left( {{y_{i_{0}}\lbrack n\rbrack} - {{cx}^{({m + 1})}\lbrack n\rbrack} - {b_{i_{0},1}^{\prime{({m + 1})}}{w_{1}^{({m + 1})}\lbrack n\rbrack}} - {b_{i_{0},2}^{\prime{({m + 1})}}{w_{2}^{({m + 1})}\lbrack n\rbrack}} - {b_{i_{0},3}^{\prime{({m + 1})}}{w_{i_{0},3}^{({m + 1})}\lbrack n\rbrack}}} \right)^{2}}}}} & (34)\end{matrix}$Iterating between equations (33) and (34) may lead to improved estimatesfor c_(i) ₀ , b′_(i) ₀ _(,1), b′_(i) ₀ _(,2), and b′_(i) ₀ _(,3).

Discrete-Time Delay Unknown

In the example above, the discrete-time delay d was assumed as known.Nevertheless, in practice it must be estimated. Let the integers d_(min)and d_(max) be the assumed minimum and maximum values for d. Then, thecomplete V-MBF algorithm follows for d=d_(min), . . . , d_(max).

1. minimize ϕ(x, d, r₂, θ_(r) ₃ , θ_(b′), c) using the algorithm in theV-MBF algorithm assuming known discrete-time delay section withinSection A (i.e., V-MBF algorithm assuming know discrete-time delay)

2. Store parameter estimates and value of least-squares objectivefunction end

The preferred estimates for the K-FADS parameters in this approachproduce the smallest least-square error. The estimates for the MBFparameters are obtained using equation (15) and the estimates:

$\begin{matrix}{\hat{\beta} = {\frac{2}{T}\frac{1 - {\hat{r}}_{2}}{1 + {\hat{r}}_{2}}}} & (35) \\{{\hat{k}}_{i} = {\frac{2}{T}\frac{1 - {\hat{r}}_{i,3}}{1 + {\hat{r}}_{i,3}}}} & (36) \\{{\hat{b}}_{i,1} = {{\hat{b}}_{i,1}^{\prime}\frac{2}{T}}} & (37) \\{{\hat{b}}_{i,2} = {{\hat{b}}_{i,2}^{\prime}\left( {\frac{2}{T} + \hat{\beta}} \right)}} & (38) \\{{{\hat{b}}_{i,3} = {{\hat{b}}_{i,3}^{\prime}\left( {\frac{2}{T} + {\hat{k}}_{i}} \right)}},} & (39)\end{matrix}$and the assumed values c_(i,2)=0.05 and c_(i,3)=0.8 for all i.

Initialization Procedure

To start the V-MBF algorithm, initial estimates for the RV factor x, r₂,and θ_(r) ₃ needed. One method for obtaining an initial estimate for theRV factor x⁽⁰⁾ is for a user selected TAC of a voxel to be used that canessentially remain in the right ventricle throughout the duration of thescan. To generate initial estimates for r₂, and θ_(r) ₃ , the knownSteiglitz-McBride algorithm may be used.

To develop a method for computing an initial estimate, firstobservations can be used from equation (22) that for some 3rd-orderpolynomial Q_(i)(z) and 2nd-order polynomial P′_(i)(z), the z-transformof the ith TAC, is given by

$\begin{matrix}{{{Y_{i}(z)} = {{c_{i}{X(z)}} + {{X(z)}\frac{{P_{i}^{\prime}(z)}\left( {1 + z^{- 1}} \right)}{Q_{i}(z)}z^{- d}}}},{i = 1},2,\ldots\mspace{11mu},I,} & (40)\end{matrix}$where the roots of Qi(z) are r₁, r₂, and r_(i,3). Alternatively, forsome polynomial P′_(i)(z), equation (40) can be written as

$\begin{matrix}{{{Y_{i}(z)} = {{X(z)}\frac{P_{i}(z)}{Q_{i}(z)}}},{i = 1},2,\ldots\mspace{11mu},I,} & (41)\end{matrix}$where the unknown numerator polynomial is of the formP_(i)(z)=p_(i,0)+p_(i,1)z⁻¹+p_(i,2)z⁻²+p_(i,3)z⁻³+p_(i,d)z^(−d)+p_(i,d+1)z^(−(d+1))+p_(i,d+2)z^(−(d+2))+p_(i,d+3)z^(−(d+3))because P_(i)(z)=c_(i)Q_(i)(z)+P′_(i)(z)(1+z⁻¹)z^(−d). In other words,each TAC is the output of an autoregressive moving-average (ARMA) modelknown in the art that is driven by the RV factor x[n].

Given an input-output pair for a linear, time-invariant system modeledas an ARMA model, the Steiglitz-McBride algorithm can provide estimatesfor ARMA parameters. Thus, given a TAC, y_(i)[n], and initial RV factorx⁽⁰⁾, the Steiglitz-McBride algorithm can be used, which is an iterativealgorithm, to estimate P_(i)(z) and Q_(i)(z). The Steiglitz-McBridealgorithm can be summarized below

$\begin{matrix}{{\left( {p_{i}^{({m + 1})},q_{i}^{({m + 1}}} \right) = {\arg\;{\min\limits_{p,q}{\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{{\frac{{Y_{i}\left( e^{j\;\omega} \right)}{Q\left( e^{j\;\omega} \right)}}{Q_{i}^{(m)}\left( e^{j\;\omega} \right)} - \frac{{X^{(0)}\left( e^{j\;\omega} \right)}{P\left( e^{j\;\omega} \right)}}{Q_{i}^{(m)}\left( e^{j\;\omega} \right)}}}^{2}d\;\omega}}}}}},} & (42)\end{matrix}$where the discrete-time Fourier transform is used to obtain:

$\begin{matrix}{{Y_{i}\left( e^{j\;\omega} \right)} = {\sum\limits_{k = 0}^{N - 1}\;{{y_{i}\lbrack k\rbrack}e^{{- j}\;\omega\; k}}}} & (43) \\{{X^{(0)}\left( e^{j\;\omega} \right)} = {\sum\limits_{k = 0}^{N - 1}\;{{x^{(0)}\lbrack k\rbrack}e^{{- j}\;\omega\; k}}}} & (44) \\{{P\left( e^{j\;\omega} \right)} = {{\sum\limits_{k = 0}^{3}\;{p_{k}e^{{- j}\;\omega\; k}}} + {\sum\limits_{k = 0}^{3}\;{p_{d + k}e^{{- j}\;\omega\; k}}}}} & (45) \\{{Q\left( e^{j\;\omega} \right)} = {1 + {\sum\limits_{k = 1}^{3}\;{q_{k}e^{{- {j\omega}}\; k}}}}} & (46) \\{p = \left\lbrack {p_{0},p_{1},p_{2},p_{3},p_{d},p_{d + 1},p_{d + 2},p_{d + 3}} \right\rbrack} & (47) \\{q = {\left\lbrack {q_{1},q_{2},q_{3}} \right\rbrack.}} & (48)\end{matrix}$

Also P_(i) ^((m))(e^(jw) ), Q_(i) ^((m))(e^(jw) ), p_(i) ^((m+1)), andq_(i) ^((m+1) can be similarly defined, and Q_(i) ⁽⁰⁾(e^(jw))=1 is thechosen initial estimate for Q(e^(jw)). Note, from Parseval's theorem,the objective function in (42) can be equivalent to a linearleast-squares objective function. Thus, at the mth iteration, theSteiglitz-McBride algorithm entails a filtering step (i.e., initial RVfactor x⁽⁰⁾[n] and ith TAC, y_(i)[n], are filtered by 1/Q_(i) ^((m))(z))and minimization of a linear least-squares objective function.

Let {circumflex over (Q)}_(i)(z) denote the resulting estimate forQ_(i)(z) using the Steiglitz-McBride algorithm andz_(i,1)≥z_(i,2)≥z_(i,3) the roots of {circumflex over (Q)}_(i)(z). Theestimates for r₂ and r_(i,3) are obtained from the roots of {circumflexover (Q)}_(i)(z), I=1, 2, . . . , I, in the following manner. Because ifβ<1 and β>k_(i) by an order of magnitude and one of the roots ofQ_(i)(z) equals one, in theory, the initial estimates for the parametersr₂ and r_(i,3) can be r₂ ⁽⁰⁾=avg{z_(1,2), z_(2,2), . . . , z_(I,2)} andr_(i,3) ⁽⁰⁾=z_(i,3), respectively.

The Steiglitz-McBride algorithm is not used instead of the V-MBFalgorithm to estimate the discrete-time K-FADS parameters when given anestimate for the RV factor because, for one reason, the roots of{circumflex over (Q)}_(i)(z) are not guaranteed to be real with one rootconstrained to equal one, as required by the discrete K-FADS model.Another reason is that simulations discovered that estimating c andθ_(b), from the estimates {circumflex over (P)}_(i)(z), {circumflex over(Q)}_(i)(z), i=1, 2, . . . , I, was not reliable.

Simulation Studies

To assess the potential of the V-MBF algorithm, simulated data wasapplied that modeled patient data used for cardiac health assessments.The model for the RV curve wasf(t)=A(t−t ₀)exp(−α(t−t ₀))u(t−t ₀),  (49)where A=2700, t₀=14 seconds, and α=0.12. Referring to equations (4),(5), (6), and (7), the parameters of the simulated data were γ=0.2,β=0.1, τ=13 seconds, and for all i, k_(i,2)=0.001 s⁻¹, k_(i,3)=0.01 s⁻¹,c_(i)=0.03, c_(i,2)=0.05, c_(i,3)=0.8, c_(i,4)=0.8. Note, the V-MBFalgorithm does not require that the parameters c_(i), k_(i,2), andk_(i,3) be voxel independent. The values for these parameters weresimply chosen for exemplary purposes. The MBF parameters {k_(i,1)}(units s⁻¹) for the 20 voxel scenario that were considered werek₁=[0.0075,0.0462,0.0371,0.0233,0.0534,0.0281,0.0320,0.0336,0.0433,0.0036,0.00830.0034,0.0021,0.0103,0.0096,0.0345,0.0316,0.0031,0.0257,0.0346].  (50)

The simulated TACs a_(i)(t), I=1, 2, . . . , I were computed usingequation (11) and the above parameter values.

To model the integration characteristic of the scanner, the simulatedTACs were integrated using a time-varying window. The resultingintegrated data simulated the activity data that would actually beavailable in practice. A typical protocol used at an exemplary PETCenter could leads to the following specification for the simulatedintegrated TACs (I-TACs)

$\begin{matrix}{g_{k} = \left\{ {\begin{matrix}{{\frac{1}{T_{1}}{\int_{{({k - 1})}T_{1}}^{{kT}_{1}}{{a_{i}(t)}\ d\; t}}},} & {{k = 1},2,\ldots\mspace{14mu},20} \\{{\frac{1}{T_{2}}{\int_{{({k - 1})}T_{2}}^{{kT}_{2}}{{a_{i}(t)}\ d\; t}}},} & {{k = 21},22,\ldots\mspace{14mu},25} \\{{\frac{1}{T_{3}}{\int_{{({k - 1})}T_{3}}^{{kT}_{3}}{{a_{i}(t)}\ d\; t}}},} & {{k = 26},27,\ldots\mspace{14mu},31}\end{matrix},} \right.} & (51)\end{matrix}$where T₁=3 seconds, T₂=12 seconds, and T₃=30 seconds.

The V-MBF algorithm can be based on standard TAC data (i.e., a_(i)(nT)).Consequently, the I-TAC data is preferrably pre-processed. The I-TACdata is assumed to be nearly piece-wise linear. It follows using a knownmethod from Kuhle that the standard TAC data at the midpoints of thewindows is approximately:a _(i)(0.5(kT ₁+(k−1)T ₁))≈g _(k) , k=1, 2, . . . , 20  (52)a _(i)(0.5(kT ₂+(k−1)T ₂))≈g _(k) , k=21, 22, . . . , 25  (53)a _(i)(0.5(kT ₃+(k−1)T ₃))≈g _(k) , k=26, 27, . . . , 31  (54)Now, standard sampled TAC data can be estimated from the measuredactivity data {g_(k)} using interpolation. Specifically, the “known”values for a_(i)(t), {a_(i)(0.5 (kT₁+(k−1)T₁))}_(k=1) ²⁰, {a_(i)(0.5(kT₂+(k−1)T₂))}_(k=21′) ²⁵, and {a_(i)(0.5 (kT₃+(k−1)T₃))}_(k=31) ²⁶ canbe used (see equations (52), (53), and (54)), and linear interpolationto obtain estimates for y_(i)[n]=a_(i)(nT), n=0, 1, . . . , N−1. In thesimulations, a preferred sampling interval is T=0.05 sec. It is notedthat the approach described above for obtaining standard sampled TACdata would also be used to generate an initial RV factor from I-TAC datalocated in the RV.

The V-MBF algorithm described above was applied for 5000 iterationswhere one sub-iteration was used to update the estimate for the RVcurve. In this simulation, the maximum error in the MBF estimates was1.5 percent. A typical result of a V-MBF algorithm is summarized in theFIG. 2, where (a) shows a true, initial, and estimated RV factor, (b)compares a true TAC (computed analytically) and with one generated usingthe discrete-time K-FADS model and estimated parameters, (c) is the sameas (b) except for a different voxel, and (d) is a plot of theleast-squares objective function as a function of iteration number. Inthe exemplary study, the V-MBF algorithm was stable and monotonicallydecreased the least-squares objective function (see FIG. 2(d)).

Thus, the V-MBF algorithm can be based on a model that accounts for thefact that the shape of TACs due to ischemic and normal tissue aredifferent. In fact, the model can allow for the factors that representfree and trapped ammonia to be voxel dependent and physiologicallyappropriate. By contrast, in a standard FADS model, it is assumed thatTACs in ischemic and normal tissue can be modeled as a linearcombination of the same three factors. The present methods and systemsrepresent a significant improvement in the art as a more appropriatemodel to provide more accurate MBF estimates than available methods.

The V-MBF algorithm presented herein performs well in simulation studieswhere unknown MBF parameters varied by an order of magnitude. Thissuggests that the V-MBF algorithm is robust and would perform well inpractice, where MBF values due to ischemic and normal tissue can varyover a wide range. Although random noise was not added to the simulatedTAC data, interpolation noise and noise due to the discretization of thecontinuous-time K-FADS model were present. Also, because only theintegrated RV factor is available, the first time point where the RVfactor is nonzero can never be known with certainty. With these threesources of noise, the maximum error of the MBF parameters estimates was1.5 percent. It should be noted that including more data points should,lead to improved MBF estimates because the parameter β is voxelindependent.

II. A Second Approach

A second approach to MBF analysis will be described as follows. Forclarity, equations will be renumbered to start with equation (1) withinthe discussion of this other approach.

As discussed above, the problem of estimating the weights and signalvectors of the above described models is a blind inverse problem. Inthis case, the data are nonnegative J×1 vectors {a₁, a₂, . . . , a_(J)}that are modeled as a weighted sum of signal vectors plus noise

$\begin{matrix}{{a_{i} = {{\sum\limits_{k = 1}^{K}\;{c_{ik}f_{k}}} + e_{i}}},\mspace{14mu}{i = 1},2,\ldots\mspace{14mu},I} & (1)\end{matrix}$where the J×1 signal vectors, {f_(k)}, and signal weights, {c_(i,k)},are nonnegative and unknown, and {e_(i)} is the noise. In our discussionof this approach, it will be convenient to use an equivalent expressionfor the model

$\begin{matrix}{{a_{ij} = {{\sum\limits_{k = 1}^{K}\;{c_{ik}f_{jk}}} + e_{ij}}},{i = 1},2,\ldots\mspace{14mu},I,{j = 0},1,2,\ldots\mspace{14mu},{J - 1},} & (2)\end{matrix}$where a_(ij) and e_(ij) are the jth components of a_(i) and e_(i),respectively, and f_(jk) is the jth value of the kth signal vector. Anexample where this model is used is cardiac imaging using dynamicpositron emission tomography (PET), where a_(ij) is a measure of theradiopharmaceutical concentration in the ith voxel at the jth timepoint. Another example is multispectral imaging, where a_(ij) representthe value of the ith pixel in the jth spectral plane.

Given the data {a_(ij)}, the problem is to estimate the weights {c_(jk)}and signal vector values {f_(jk)}. The least-squares estimates of theweights and signal vector values are obtained by minimizing theleast-squares objective function L subject to a non-negativityconstraint.

$\begin{matrix}{{\left( {\hat{c},\hat{f}} \right) = {\arg{\min\limits_{{c \geq 0},{f \geq 0}}{L\left( {c,f} \right)}}}},{where}} & (3) \\{{L\left( {c,f} \right)} = {\sum\limits_{i = 1}^{I}\;{\sum\limits_{j = 0}^{J - 1}\;{\left( {a_{ij} - {\sum\limits_{k = 1}^{K}\;{c_{ik}f_{jk}}}} \right)^{2}.}}}} & (4)\end{matrix}$The vectors c and f contain the signal weights and signal vector values,respectively:c

[c₁₁, c₂₁, . . . , c_(I1) , c ₁₂, c₂₂, . . . , c_(I2), . . . , c_(1K),c_(2K), . . . , c_(IK)]  (5)f

[f₁₁, f₂₁, . . . , f_(I1), f₁₂, f₂₂, . . . , f_(I2), . . . , f_(1K),f_(2K), . . . , f_(IK)]  (6)The least squares estimation problem is ill-posed, so the results arehighly dependent on the initial estimates.

Estimation of Signal Vectors and Weights

A known standard least-squares algorithm by Lee and Seung monotonicallydecreased the least-squares objective function, L, and producednonnegative estimates. In this section, this “standard least-squaresalgorithm” is re-derived using a known technique called themajorize-minimize method (MM). This derivation will place the proposedleast-squares extensions in context so that their advantages over thestandard least-squares algorithm are clear.

An approach for minimizing the least-squares objective function, L,would be, in an alternating fashion, to minimize L with respect to thesignal weights while holding the signal vectors fixed to their currentvalue, and then minimize L with respect to the signal vectors whileholding the signal weights fixed to their current value. Given initialestimates {c_(ik) ⁽⁰⁾} and {f_(jk) ⁽⁰⁾}, this algorithm can be expressedmathematically as follows:

Step 1′$c^{({n + 1})} = {\arg\mspace{11mu}{\min\limits_{c \geq 0}\;{L\left( {c,f^{(n)}} \right)}}}$Step 2′$f^{({n + 1})} = {\arg\mspace{11mu}{\min\limits_{f \geq 0}\;{L\left( {c^{({n + 1})},f} \right)}}}$Step 3′ Repeat Steps 1′ and 2′ for a fixed number of iterations or untilsome desired stopping criterion is metHere, Steps 1′ and 2′ imply that the resulting algorithm monotonicallydecreases the least-squares objective functionL(c ^((n+1)) ,f ^((n+1)))≤L(c ^((n)) ,f ^((n))) for all n=0,1,2, . . .  (7)

Introduction to the MM Method

The minimizations problems in Steps 1′ and 2′ are difficult so analternative approach is needed. In this section, the MM method isintroduced: then, the following section demonstrates how the MM methodcan be used to develop an algorithm that, by construction, monotonicallydecreases the least-squares objective function, produces nonnegativeestimates, and is straightforward to implement.

Consider a real valued function f with domain D∈R^(n) that is to beminimized. A real valued function g with domain {(x₁, x₂): x₁, x₂∈D} issaid to majorize the function f if the following conditions hold for allx, y∈D:g(x,y)≥f(x) and  (C1′)g(x,x)=f(x).  (C2′)MM algorithms are a viable approach provided a majorizing function canbe found that is easier to minimize than the original objectivefunction. Assuming that g is a suitable majorizing function for f, thecorresponding MM algorithm is

$\begin{matrix}{x^{({k + 1})} = {\arg{\min\limits_{x \in D}{g\left( {x,x^{(k)}} \right)}}}} & (8)\end{matrix}$It follows from (C1′) and (C2′) that the MM algorithm defined by (8) ismonotonicf(x ^((k+1)))≤g(x ^((k+1)) ,x ^((k)))≤g(x ^((k)) ,x ^((k)))=f(x^((k))).  (9)

Standard Least-Squares Algorithm

The following describes an MM algorithm for estimating signal vectorsand signal weights that is equivalent to the least-squares algorithm ofLee and Seung.

In the discussion that follows, it will be convenient to define setsD_(c) and D_(f), where D_(c) and D_(f) are the set of non-negativevectors of dimension IK×1 and JK×1, respectively. Let c^((n)) andf^((n)) be the current estimates for the signal weights and signalvectors, respectively. Further, let q and r be certain majorizingfunctions that satisfy the following conditions for x, y, c∈D_(c) and s,t, f∈D_(f) beq(x,y,f)≥L(x,f)  (C1)q(x,x,f)=L(x,f)  (C2)r(s,t,c)≥L(c,s)  (C3)r(x,x,c)=L(c,x).  (C4)Using the idea behind equation (8), CLA put forth the following MMalgorithm for minimizing the least-squares objective function L

$\begin{matrix}{c^{({n + 1})} = {\arg{\min\limits_{c \geq 0}{q\left( {c,c^{(n)},f^{(n)}} \right)}}}} & (10) \\{f^{({n + 1})} = {\arg{\min\limits_{f \geq 0}{{r\left( {f,f^{(n)},c^{({n + 1})}} \right)}.}}}} & (11)\end{matrix}$It is straightforward to show that least-squares objection functionmonotonically decreases with increasing iterations. First, from (C1) andequation (10) it follows, respectively thatL(c ^((n+1)) ,f ^((n)))≤q(c ^((n+1)) ,c ^((n)) ,f ^((n)))  (12)q(c ^((n+1)) ,c ^((n)) ,f ^((n)))≤q(c ^((n)) ,c ^((n)) ,f ^((n))).  (13)Similarly, from (C3) and equation (11) it follows, respectively, thatL(c ^((n+1)) ,f ^((n+1)))≤r(f ^((n+1)) ,f ^((n)) ,c ^((n+1)))  (14)r(f ^((n+1)) ,f ^((n)) ,c ^((n+1)))≤r(f ^((n)) ,f ^((n)) ,c^((n+1))).  (15)Now, from (C2) and (C4) it follows, respectively, thatq(c ^((n)) ,c ^((n)) ,f ^((n)))=L(c ^((n)) ,f ^((n)))  (16)r(f ^((n)) ,f ^((n)) ,c ^((n+1)))=L(c ^((n+1)) ,f ^((n))).  (17)Consequently, we can conclude from equations (12)-(17) thatL(c ^((n+1)) ,f ^((n+1)))≤L(c ^((n)) ,f ^((n))).  (18)

At this point, all that remains is to determine the majorizing functionsq and r. From equation (4), the least-squares objection function can bewritten as

$\begin{matrix}{{L\left( {c,f} \right)} = {\sum\limits_{i = 1}^{L}\;{\sum\limits_{j = 0}^{J - 1}\;{\left\lbrack {a_{ij}^{2} - {2\; a_{ij}{\sum\limits_{k = 1}^{K}\;{c_{ik}f_{jk}}}} + \left( {\sum\limits_{k = 1}^{K}\;{c_{ik}f_{jk}}} \right)^{2}} \right\rbrack.}}}} & (19)\end{matrix}$Using a known approach, the convexity of the square function isexploited to obtain the following inequality:

$\begin{matrix}{\left( {\sum\limits_{k = 1}^{K}\;{c_{ik}f_{jk}}} \right)^{2} = \left( {\sum\limits_{k = 1}^{K}\;{\frac{c_{ik}^{(n)}f_{jk}}{\left\lbrack {\sum\limits_{k^{\prime} = 1}^{K^{\prime}}\;{c_{{ik}^{\prime}}^{(n)}f_{{jk}^{\prime}}}} \right\rbrack}\frac{\left\lbrack {\left\lbrack {\sum\limits_{k^{\prime} = 1}^{K^{\prime}}\;{c_{{ik}^{\prime}}^{(n)}f_{{jk}^{\prime}}}} \right\rbrack c_{ik}} \right\rbrack}{c_{ik}^{(n)}}}} \right)^{2}} & (20) \\{\leq {\sum\limits_{k = 1}^{K}\;{\frac{c_{ik}^{(n)}f_{jk}}{\left\lbrack {\sum\limits_{k^{\prime} = 1}^{K^{\prime}}\;{c_{{ik}^{\prime}}^{(n)}f_{{jk}^{\prime}}}} \right\rbrack}\left( \frac{\left\lbrack {\sum\limits_{k^{\prime} = 1}^{K^{\prime}}\;{c_{{ik}^{\prime}}^{(n)}f_{{jk}^{\prime}}}} \right\rbrack c_{ik}}{c_{ik}^{(n)}} \right)^{2}}}} & (21) \\{{\leq {{{\hat{a}}_{ij}\left( {c^{(n)},f} \right)}{\sum\limits_{k = 1}^{K}\;\left( {\frac{c_{ik}^{2}}{c_{ik}^{(n)}}f_{jk}} \right)}}},} & (22)\end{matrix}$where c^((n))∈D_(c), f∈D_(f), and

$\begin{matrix}{{{\hat{a}}_{ij}\left( {c,f} \right)}\overset{\Delta}{=}{\sum\limits_{k = 1}^{K}\;{c_{ik}{f_{jk}.}}}} & (23)\end{matrix}$It should be noted that (20) is a convex combination with weightsw_(ijk)

c_(ik) ^((n)) f_(jk)[Σ_(k′=1) ^(K′)c_(ik′) ^((n))f_(jk′)]⁻¹, wherew_(ijk)≥0 and Σ_(k=1) ^(K) w_(ijk)=1 for all i, j. Replacing the(Σ_(k=1) ^(K) c_(ik)f_(jk))² term in (19) by the right hand side of(22), we get the following majorizing function for L

$\begin{matrix}{{{q\left( {c,c^{(n)},f} \right)}\overset{\Delta}{=}{\sum\limits_{i = 1}^{I}\;{\sum\limits_{j = 0}^{J - 1}\;\left\{ {a_{ij}^{2} - {2\; a_{ij}{\sum\limits_{k = 1}^{K}\;{c_{ik}f_{jk}}}} + {{{\hat{a}}_{ij}\left( {c^{(n)},f} \right)}{\sum\limits_{k = 1}^{K}\;\left( {\frac{c_{ik}^{2}}{c_{ik}^{(n)}}f_{jk}} \right)}}} \right\}}}},} & (24)\end{matrix}$where, by construction, q(c, c^((n)), f)≥L(c, f) and q(c, c, f)=L(c, f)for all c, c^((n))∈D_(c) and f∈D_(f).

By repeating the steps used to derive equation (24) with the roles of cand f switched, the majorizing function

$\begin{matrix}{{{r\left( {f,f^{(n)},c} \right)}\overset{\Delta}{=}{\sum\limits_{i = 1}^{I}\;{\sum\limits_{j = 0}^{J - 1}\;\left\{ {a_{ij}^{2} - {2\; a_{ij}{\sum\limits_{k = 1}^{K}\;{c_{ik}f_{jk}}}} + {{{\hat{a}}_{ij}\left( {c,f^{(n)}} \right)}{\sum\limits_{k = 1}^{K}\;\left( {\frac{f_{ik}^{2}}{f_{ik}^{(n)}}c_{jk}} \right)}}} \right\}}}},} & (25)\end{matrix}$which satisfies the properties r(f, f^((n)), c)≥L(c, f) and r(f, f,c)=L(c, f) for all c∈D_(c) and f, f^((n))∈D_(f).

To determine updates defined in equations (10) and (11), the partialderivatives of q(c, c^((n)), f^((n))) and r(f, f^((n)), c^((n+1))) arecomputed with respect to c and f, respectively, with the correspondingequations set to zero. It is straightforward to show that thederivatives are given by

$\begin{matrix}{\mspace{79mu}{\frac{\partial{q\left( {c,c^{(n)},f^{(n)}} \right)}}{\partial c_{i^{\prime}k^{\prime}}} = {\sum\limits_{j = 0}^{J - 1}\;\left\lbrack {{{- 2}\; a_{i^{\prime}j}f_{{jk}^{\prime}}^{(n)}} + {2\frac{c_{i^{\prime}k^{\prime}}}{c_{i^{\prime}k^{\prime}}^{(n)}}f_{{jk}^{\prime}}^{(n)}{{\hat{a}}_{ij}\left( {c^{(n)},f^{(n)}} \right)}}} \right\rbrack}}} & (26) \\{\frac{\partial{r\left( {f,f^{(n)},c^{({n + 1})}} \right)}}{\partial f_{j^{\prime}k^{\prime}}} = {\sum\limits_{i = 1}^{I}\;{\left\lbrack {{{- 2}\; a_{i^{\prime}j^{\prime}}c_{{ik}^{\prime}}^{({n + 1})}} + {2\frac{f_{j^{\prime}k^{\prime}}}{f_{j^{\prime}k^{\prime}}^{(n)}}c_{{ik}^{\prime}}^{({n + 1})}{{\hat{a}}_{ij}\left( {c^{({n + 1})},f^{(n)}} \right)}}} \right\rbrack.}}} & (27)\end{matrix}$Setting the derivatives in equations (26) and (27) to zero leads to thefollowing least-squares algorithm for estimating the signal weights andsignal vectors:

$\begin{matrix}{{c_{ik}^{({n + 1})} = {c_{ik}^{(n)}\frac{\sum\limits_{j = 0}^{J - 1}\;{f_{jk}^{(n)}a_{ij}}}{\sum\limits_{j = 0}^{J - 1}\;{f_{jk}^{(n)}{{\hat{a}}_{ij}\left( {c^{(n)},f^{(n)}} \right)}}}}},{i = 1},2,\ldots\mspace{14mu},I,{k = 1},2,\ldots\mspace{14mu},K} & (28) \\{{f_{jk}^{({n + 1})} = {f_{jk}^{(n)}\frac{\sum\limits_{i = 0}^{I - 1}\;{c_{ik}^{({n + 1})}a_{ij}}}{\sum\limits_{i = 0}^{I - 1}\;{c_{ik}^{(n)}{{\hat{a}}_{ij}\left( {c^{({n + 1})},f^{(n)}} \right)}}}}},{j = 0},1,2,\ldots\mspace{14mu},I,{k = 1},2,\ldots\mspace{14mu},{K.}} & (29)\end{matrix}$

As shown above, the starting point for developing the proposed MMalgorithm was equations (10) and (11). Alternatively, an MM algorithmcan be developed by updating the signal vectors first and then thesignal weights:

$\begin{matrix}{f^{({n + 1})} = {\arg{\min\limits_{f \geq 0}{r\left( {f,f^{(n)},c^{(n)}} \right)}}}} & (30) \\{c^{({n + 1})} = {\arg{\min\limits_{c \geq 0}{{q\left( {c,c^{(n)},f^{({n + 1})}} \right)}.}}}} & (31)\end{matrix}$In this case, the resulting MM algorithm is

$\begin{matrix}{{f_{jk}^{({n + 1})} = {f_{jk}^{(n)}\frac{\sum\limits_{i = 1}^{I}{c_{ik}^{(n)}a_{ij}}}{\sum\limits_{i = 1}^{I}\;{c_{ik}^{(n)}{{\hat{a}}_{ij}\left( {c^{(n)},f^{(n)}} \right)}}}}},{j = 0}, 1, 2,\ldots\mspace{14mu},{J - 1},{k = 1}, 2,\ldots\mspace{14mu},K} & (32) \\{{c_{ik}^{({n + 1})} = {c_{ik}^{(n)}\frac{\sum\limits_{j = 0}^{J - 1}{f_{jk}^{({n + 1})}a_{ij}}}{\sum\limits_{j = 0}^{J - 1}{f_{jk}^{({n + 1})}{{\hat{a}}_{ij}\left( {c^{(n)},f^{({n + 1})}} \right)}}}}},{i = 1}, 2,\ldots\mspace{14mu}, I,{k = 1}, 2,\ldots\mspace{14mu},{K.}} & (33)\end{matrix}$

Additive Least-Squares Algorithm

The standard least-squares algorithm defined by equations (28) and (29)(see also (32) and (33)) is referred to as a multiplicative algorithm.In this section, we derive an additive algorithm by developing analternative majorizing function for the least-squares objective functionL.

Exploiting the convexity of the square function we obtain the followinginequality

$\begin{matrix}{\left( {\sum\limits_{k = 1}^{K}\;{c_{ik}f_{jk}}} \right)^{2} = \left\lbrack {\sum\limits_{k = 1}^{K}\;{\frac{f_{jk}}{{\overset{\_}{f}}_{j}}\left( {{{\overset{\_}{f}}_{j}c_{ik}} - {{\overset{\_}{f}}_{j}c_{ik}^{(n)}} + {{\hat{a}}_{ij}\left( {c^{(n)},f} \right)}} \right)}} \right\rbrack^{2}} & (34) \\{\leq {\sum\limits_{k = 1}^{K}\;{\frac{f_{jk}}{{\overset{\_}{f}}_{j}}\left( {{{\overset{\_}{f}}_{j}c_{ik}} - {{\overset{\_}{f}}_{j}c_{ik}^{(n)}} + {{\hat{a}}_{ij}\left( {c^{(n)},f} \right)}} \right)^{2}}}} & (35)\end{matrix}$where f _(j)

Σ_(k=1) ^(K) f_(jk), c^((n))∈D_(c), and f∈D_(f). It should be noted that(34) is a convex combination with weights w_(jk)

f_(jk)/f _(j), where w_(jk)≥0 for all j,k and Σ_(k=1) ^(K) w_(jk)=1 forall j. Replacing the (Σ_(k=1) ^(K) c_(ik)f_(jk))² term in (19) by theright hand side of (35), we get the following majorizing function for L

$\begin{matrix}{{{q_{A}\left( {c,c^{(n)},f} \right)}\overset{\Delta}{=}{\sum\limits_{i = 1}^{I}\;{\sum\limits_{j = 0}^{J - 1}\;\left\{ {a_{ij}^{2} - {2\; a_{ij}{\sum\limits_{k = 1}^{K}\;{c_{ik}f_{jk}}}} + {\sum\limits_{k = 1}^{K}\;{\frac{f_{jk}^{2}}{f_{jk}^{(n)}}\left( {{{\overset{\_}{f}}_{j}c_{jk}} - {{\overset{\_}{f}}_{j}c_{ik}^{(n)}} + {{\hat{a}}_{ij}\left( {c^{(n)},f} \right)}} \right)^{2}}}} \right\}}}},} & (36)\end{matrix}$where, by construction, q_(A)(c, c^((n)), f)≥L(c, f) and q_(A)(c, c,f)=L(c, f) for all c, c^((n))∈D_(c) and f∈D_(f).

When the steps used to derive equation (24) are repeated with the rolesof c and f switched, we get the following majorizing function:

$\begin{matrix}{{{r_{A}\left( {f,f^{(n)},c} \right)}\overset{\Delta}{=}{\sum\limits_{i = 1}^{I}\;{\sum\limits_{j = 0}^{J - 1}\;\left\{ {a_{ij}^{2} - {2a_{ij}{\sum\limits_{k = 1}^{K}\;{c_{ik}f_{jk}}}} + {\sum\limits_{k = 1}^{K}\;{\frac{c_{ik}}{{\overset{\_}{c}}_{i}}\left( {{{\overset{\_}{c}}_{i}f_{jk}} - {{\overset{\_}{c}}_{i}f_{jk}^{(n)}} + {{\hat{a}}_{ij}\left( {c,f^{(n)}} \right)}} \right)^{2}}}} \right\}}}},} & (37)\end{matrix}$which satisfies the properties, r_(A)(f, f^((n)), c)≥L(c, f) andr_(A)(f, f, c)=L (c, f) for all c∈D_(c) and f, f^((n))∈D_(f).

The updates defined by equations (10) and (11) are determined with q andr replaced by q_(A) and r_(A), respectively. Specifically, the partialderivatives of q_(A)(c, c^((n)), f^((n))) and r_(A)(f, f^((n)),c^((n+1))) are computed with respect to c and f, respectively, and thenthe resulting equations are set to zero. The desired derivatives can begiven by

$\begin{matrix}{{\frac{\partial{q_{A}\left( {c,c^{(n)},f^{(n)}} \right)}}{{\partial c_{i}},_{k},} = {{- 2}{\sum\limits_{j = 0}^{J - 1}\; a_{i}}}},_{j}{f_{{jk}^{\prime}}^{(n)} + {2{\sum\limits_{j = 0}^{J - 1}\;{{\frac{f_{{jk}^{\prime}}^{(n)}}{{\overset{\_}{f}}_{j}^{(n)}}\left\lbrack {{{\overset{\_}{f}}_{j}^{(n)}c_{i}},_{k},{{{- {\overset{\_}{f}}_{j}^{(n)}}c_{i^{\prime}k^{\prime}}^{(n)}} + {\hat{a}}_{i}},_{j}\left( {c^{(n)},f^{(n)}} \right)} \right\rbrack}{\overset{\_}{f}}_{j}^{(n)}}}}}} & (38) \\{{\frac{\partial{r_{A}\left( {f,f^{(n)},c^{({n + 1})}} \right)}}{{\partial f_{i}},_{k},} = {{- 2}{\sum\limits_{i = 0}^{I}\; a_{i}}}},_{j}{f_{{jk}^{\prime}}^{({n + 1})} + {2{\sum\limits_{i = 0}^{i - 1}\;{{\frac{c_{{ik}^{\prime}}^{({n + 1})}}{{\overset{\_}{c}}_{i}^{({n + 1})}}\left\lbrack {{{\overset{\_}{c}}_{i}^{({n + 1})}f_{j}},_{k},{{{- {\overset{\_}{c}}_{i}^{({n + 1})}}f_{j^{\prime}k^{\prime}}^{(n)}} + {\hat{a}}_{ij}},\left( {c^{({n + 1})},f^{(n)}} \right)} \right\rbrack}{\overset{\_}{c}}_{i}^{({n + 1})}}}}}} & (39)\end{matrix}$Setting the derivatives in equations (38) and (39) to zero leads to thefollowing additive least-squares algorithm for estimating the signalweights and signal vectors

$\begin{matrix}{{c_{ik}^{({n + 1})} = \left\lbrack {c_{ik}^{(n)} + \frac{\sum_{j = 0}^{J - 1}{f_{jk}^{(n)}\left\lbrack {a_{ij} - {{\hat{a}}_{ij}\left( {c^{(n)},f^{(n)}} \right)}} \right\rbrack}}{\sum_{j = 0}^{J - 1}{f_{jk}^{(n)}{\overset{\_}{f}}_{j}^{(n)}}}} \right\rbrack_{+}},{i = 1}, 2,\ldots\mspace{14mu}, l,{k = 1},2,\ldots\mspace{14mu},K} & (40) \\{{f_{ik}^{({n + 1})} = \left\lbrack {f_{ik}^{(n)} + \frac{\sum_{i = 1}^{i}{c_{jk}^{({n + 1})}\left\lbrack {a_{ij} - {{\hat{a}}_{ij}\left( {c^{({n + 1})},f^{(n)}} \right)}} \right\rbrack}}{\sum_{j = 0}^{J - 1}{c_{ik}^{({n + 1})}{\overset{\_}{c}}_{i}^{({n + 1})}}}} \right\rbrack_{+}},{j = 0}, 1,\ldots\mspace{14mu},{J - 1},{k = 1},2,\ldots\mspace{14mu},K,} & (41)\end{matrix}$where [x]+

max(0, x). Now, when the signal vectors are updated first and then thesignal weights, the corresponding additive least-squares algorithm isgiven by

$\begin{matrix}{{f_{jk}^{({n + 1})} = \left\lbrack {f_{jk}^{(n)} + \frac{\sum_{i = 1}^{J}{c_{ik}^{(n)}\left\lbrack {a_{ij} - {{\hat{a}}_{ij}\left( {c^{(n)},f^{(n)}} \right)}} \right\rbrack}}{\sum_{i = 1}^{I}{c_{ik}^{(n)}{\overset{\_}{c}}_{i}^{(n)}}}} \right\rbrack_{+}},{j = {0{\quad{, 1,\ldots\mspace{14mu},{J - 1},{k = 1},2,\ldots\mspace{14mu},K,}}}}} & (42) \\{{c_{ik}^{({n + 1})} = \left\lbrack {c_{ik}^{(n)} + \frac{\sum_{j = 0}^{J - 1}{f_{jk}^{({n + 1})}\left\lbrack {a_{ij} - {{\hat{a}}_{ij}\left( {c^{(n)},f^{({n + 1})}} \right)}} \right\rbrack}}{\sum_{j = 0}^{J - 1}{f_{jk}^{({n + 1})}{\overset{\_}{f}}_{j}^{({n + 1})}}}} \right\rbrack_{+}},{i = 1}, 2,\ldots\mspace{14mu}, I,{k = 1},2,\ldots\mspace{14mu},{K.}} & (43)\end{matrix}$

Extensions of Standard Least-Squares Algorithm: Application toMyocardial Blood Flow Estimation Using Positron Emission Tomography

In this section, the standard least-squares algorithm application to theproblem of estimating absolute myocardial blow (MBF) noninvasively usingpositron emission tomography is addressed. Then, extensions of thestandard least-squares algorithm are presented that lead to improved MBFestimates. These extensions are also applicable to the additiveleast-squares algorithm.

To assess the heart of a patient, it is desirable to estimate thepatient's MBF noninvasively. One way to obtain this information is tofirst perform a dynamic PET scan of the patient's heart and then applyestimation algorithms that are based on equation (2) and other availablemodels for dynamic cardiac PET data. Note, in the PET literature, theterms factor curves and factor weights are used for the terms signalvectors and signal weights. For PET based MBF estimation, theleast-squares algorithm in equations (28) and (29) could first be usedto estimate the factor curves and weights for a given dynamic PET dataset. Then, using the resulting estimates, standard methods could be usedto estimate the absolute myocardial blood flow of the patient. Theaccuracy of the MBF estimates would greatly depend on the performance ofthe standard least-squares algorithm, which, as mentioned previously, ishighly dependent on the initial estimates. Therefore, we developextensions of the least-squares method that greatly reduce the parameterspace by incorporating a priori information. In practice, the proposedalgorithms are expected to be more stable and produce more accurateestimates of the factor curves and weights than the standardleast-squares algorithm. Therefore, improved MBF estimation isanticipated when the proposed algorithms are used instead of thestandard least-squares algorithm.

Model

Let a_(i)(t) denote the continuous-time activity in voxel i at time t,where t∈[0, T] and T is the duration of the scan. In practice, onlysamples of the data are available so a_(ij) denotes the activity invoxel i at time t=jT_(s)(i.e., time frame j)a _(ij)

a _(i)(jT _(s)), i=1,2, . . . , I,j=0,1,2, . . . , J−1  (44)where T_(s) is the sampling interval, I is the number of voxels, and Jis the number of time frames. Referring to equation (2), it is typicallyassumed that the activity a_(ij) is a linear combination of K=3 unknownfactor curves representing the sampled right ventricular blood pool,left ventricular blood pool, and myocardial tissue curves, respectively.Note, the factor weights for a particular time frame can be viewed as animage and therefore are collectively referred to as a factor image.

Physiological Based Constraints: Right and Left Ventricle Tissue Curvesgo to Zero

Due to the physiology of the heart and half-life of theradiopharmaceuticals used in nuclear cardiology, the factor curves forthe right and left ventricles go to zero as t approaches 0, which apriori information is incorporated into the estimation problem. Thus, analternative to the least-squares formulation is to add a penalty term tothe least-squares objective function that “forces” the factor curves forthe right ventricle f₁

[f₁₁, f₂₁, . . . , f_(J1)] and left ventricle f₂

[f₁₂, f₂₂, . . . , f_(J2)] to go to zero. Given its potential advantageover the standard least-squares method, we propose the followingpenalized least-squares method

$\begin{matrix}{{\left( {\overset{\_}{c},\hat{f}} \right) = {{\arg\mspace{14mu}{\min\limits_{{c \geq 0},{f \geq 0}}\mspace{14mu}{L\left( {c,f} \right)}}} + {\beta_{1}\Lambda\;\left( f_{1} \right)} + {\beta_{2}{\Lambda\left( f_{2} \right)}}}},} & (45)\end{matrix}$where the penalty parameters β₁ and β₂ control the degree of influenceof the penalty terms Λ(f₁) and Λ(f₂), respectively. Although there aremany possible choices for the penalty function Λ, in this approach thefollowing function is used:Λ(g)=(Σ_(j>j) ₀ ^(J−1) g _(j))²  (46)where g is a real J×1 vector. It can be seen that Λ(g) is the energy ofg after the user chosen time frame j₀.

The desired MM algorithm is obtained by setting the partial derivativesof r(f, f^((n)), c^((n+1)))+β₁Λ(f₁)+β₂ Λ(f₂) with respect to f to zero.For k=1, 2, the partial derivatives of the penalty function Λ are

$\begin{matrix}{\frac{\partial{\Lambda\left( f_{k} \right)}}{{\partial f_{j}},_{k},} = \left\{ \begin{matrix}{{2\; f_{j}},_{k},} & {{j^{\prime} > j_{0}},{k^{\prime} = k}} \\{0,} & {{otherwise},}\end{matrix} \right.} & (47)\end{matrix}$Using this result and equation (27), the following iterative algorithmis derived for estimating the right and left ventricle tissue curves

$\begin{matrix}{{f_{jk}^{({n + 1})} = {f_{jk}^{(n)}\frac{\sum_{i = 1}^{J}{c_{ik}^{({n + 1})}a_{ij}}}{{\sum_{i = 1}^{I}{c_{ik}^{({n + 1})}{{\hat{a}}_{ij}\left( {c^{({n + 1})},f^{(n)}} \right)}}} + {\beta_{k}f_{jk}^{(n)}{I_{S}(j)}}}}},{j = 0}, 1, 2,\ldots\mspace{14mu},{J - 1},{k = 1},2} & (48)\end{matrix}$where I is the indicator function and S={jo+1, jo+2, . . . , J} (i.e.,I_(s)(j)=1 for j∈S and I_(s)(j)=0 for j

S). The updates for the factor weights and tissue curve (i.e., f₃) aregiven by equations (28) and (29), respectively. The proposed algorithm,refer to herein as the PLS algorithm, monotonically decreases thepenalized least-squares (PLS) objective function in equation (45) and isguaranteed to produce nonnegative estimates for the values of the factorcurves, f, and factor weights, c. It should be noted that the PLSalgorithm provides least-squares estimates when β₁=β₂=0. Also, theupdate for the factor weights is the same for both the LS and PLSalgorithms.

Due to the nonnegativity of the factor curves for the right and leftventricles, another penalty function optionally can be used to accountfor the fact that they decrease for t sufficiently large isΛ(g)=(Σ_(j>j) ₀ ^(J−1) g _(j))²  (49)Also, with a small modification, the PLS framework can incorporate otherknown penalty functions.

Exploit Fact that Maximum of Right and Left Ventricle Tissue Curves canbe Reliably Estimated

When the radiopharmaceutical is administered to a patient, thephysiology of the body is such that the activity shows up in the rightventricle first, then the left ventricle, and finally the myocardium.These delays are such that the maximum values of the right and leftventricles are essentially free of activity from the myocardium, despitethe motion of the heart and point spread function of the PET system.Thus, an estimate of the maximum value of the right ventricle can beestimated by averaging the maximum values of TACs that lie in thecentral region of the right ventricle. In a similar way, an estimate ofthe maximum value of the left ventricle can be obtained. Methods areavailable that identify the voxels that lie in the right and leftventricles, and myocardium.

Let μ₁ and μ₂ represent the unknown maximum values of the right and leftventricles, respectively, and j₁ and j₂ denote the locations of themaximum values of the right and left ventricles. To incorporateknowledge of μ₁ and μ₂, we estimate the tissue factors usingf ^((n+1))=arg min r(f,f ^((n)) ,c ^((n))) subject to f≥0, |f _(j) ₁₁−{circumflex over (μ)}₁ |≤∈, |f _(j) ₂ ₂−{circumflex over(μ)}₂|≤∈,  (50)where ∈ is a tolerance parameter chosen by the user, and {circumflexover (μ)}₁ and {circumflex over (μ)}₂ are estimates of the maximumvalues μ₁ and μ₂.

The function r(f, f^((n)), c^((n))) is decoupled in the sense that thereare no terms of the form f_(jk)f_(j′k′), except when j=j′ and k=k′.Thus, the solution to the optimization problem in equation (50) isstraightforward and leads to the following update for the right and leftventricle tissue factors (i.e., k=1, 2):

$\begin{matrix}{{f_{jk}^{({n + 1})} = {f_{jk}^{(n)}\frac{\sum\limits_{i = 1}^{I}\;{c_{ik}^{({n + 1})}a_{ij}}}{\sum\limits_{i = 1}^{I}\;{c_{ik}^{({n + 1})}{{\hat{a}}_{ij}\left( {c^{({n + 1})},f^{(n)}} \right)}}}}},{{{for}\mspace{14mu}{all}\mspace{14mu} j} \neq j_{k}},{k = 1},2} & (51) \\{{f_{jk}^{({n + 1})} = \left\lbrack {f_{jk}^{(n)}\frac{\sum\limits_{i = 1}^{I}\;{c_{ik}^{({n + 1})}a_{ij}}}{\sum\limits_{i = 1}^{I}\;{c_{ik}^{({n + 1})}{{\hat{a}}_{ij}\left( {c^{({n + 1})},f^{(n)}} \right)}}}} \right\rbrack_{\mu_{k} - c}^{\mu_{k} + c}},{{{for}\mspace{14mu} j} = {\quad{{jk},{k = 1},{2{where}}}}}} & (52) \\{\lbrack a\rbrack_{l}^{\mu}\overset{\Delta}{=}\left\{ \begin{matrix}{a,} & {l \leq a \leq u} \\{l,} & {a < l} \\{u,} & {a > u}\end{matrix} \right.} & (53)\end{matrix}$The updates for the factor weights and tissue curve are given byequations (28) and (29), respectively. Also, the denominator terms inequations (51) and (52) would be Σ_(i=1) ^(I)c_(ik)^((n+1))â_(ij)(c^((n+1)), f^((n)))+β_(k)f_(jk) ^((n))I_(S)(j) if weincluded the penalties on the right and left ventricle tissue curvesthat were discussed in the previous section.

Reduce Unknown Parameters Via a Suitable Model for Left Ventricle TissueCurve

It has been postulated that the left ventricle tissue curve can bemodeled as the convolution of the right ventricle tissue curve with agamma function. Described below is an exploitation of this idea anddevelopment of an extension of the standard least-squares algorithm thathas significantly fewer unknowns.

1) Model: Let r(t) and l(t) denote the unknown continuous-time tissuecurves for the right and left ventricles, respectively. The leftventricle tissue curve is modeled as the convolution of the right tissuecurve with a delayed, gamma function h _(c)(t)l(t)=r(t)* h _(c)(t−τ),  (54)where “*” denotes the convolution operator

${{{\overset{\_}{h}}_{c}(t)} = {a\;\frac{t^{m - 1}}{\left( {m - 1} \right)l}{\exp\left( {- {bt}} \right)}{u(t)}\mspace{14mu}\left( {{a > 0},{b > 0},{{{and}\mspace{14mu} m} = 1},2,\ldots}\mspace{14mu} \right)}},$u(t) is the unit step function, and τ>0 is the delay. It should be notedthat the delay is due to the fact that the radiopharmaceutical activityfirst appears in the right ventricle and then the left ventricle. From(54) it follows that L(s)=R(s)H _(c)(s)exp(−sτ), where

$\begin{matrix}{{{\overset{\_}{H}}_{c}(s)}\overset{\Delta}{=}{\frac{a}{\left( {s + b} \right)^{m}}.}} & (55)\end{matrix}$

As noted previously, we let f_(j1) and F_(j2), j=0, 1, 2, . . . , J−1denoted the sampled right and left ventricle tissue curves. Where thescan durations for dynamic sequence PET protocols are not uniform, theassumption that the TACs are sampled uniformly is inappropriate.However, uniform samples of the TACs can be obtained from non-uniformsamples of the TACs via a suitable interpolation. It is of interest todetermine a discrete-time system with the property that its response tothe right ventricle tissue curve f₁ closely approximates the leftventricle tissue curve f₂. The bilinear transformation is a popular wayto transform a linear time-invariant continuous-time system into alinear time-invariant discrete-time system. A limitation of the bilineartransform is that a delay in a continuous-time system must be an integermultiple of the sampling interval. Assuming the delay τ is a multiple ofthe sampling interval Ts, the system function of the desired discretetime system, H(z), can be obtained by applying the bilineartransformation to the continuous-time system H_(c)(s)=H _(c)(s)exp(−sτ)³

$\begin{matrix}{{H(z)} = {{H_{c}(s)}❘_{s = {\frac{2}{T}\frac{1 - z^{- 2}}{1 + z^{- 1}}}}}} & (56) \\{= {{a\left\lbrack {\frac{1}{\left( {s + b} \right)^{m}}❘_{s = {\frac{2}{T}\frac{1 - z^{- 2}}{1 + z^{- 1}}}}} \right\rbrack}z^{- d}}} & (57) \\{{= {{g\left( \frac{1 + z^{- 1}}{1 - {pz}^{- 1}} \right)}^{m}z^{- d}}},} & (58)\end{matrix}$where τ=dT_(s) for some inter d, g=a(2/T_(s)+b)^(−(m+1)), andp=(2/T_(s)+b)⁻¹(2/T_(s)−b) (note: h_(c)(t)

h _(c)(t−τ)).

Let hj(θ) denote the inverse z-transform of H(z), using a notation thatillustrates its dependence on the parameters θ=(g, p, m, d). Then, theassumed relationship between the sampled right ventricle and leftventricle tissue curves can be written asf _(j2) =f _(j1) *h _(j)(θ)=Σ_(s=0) ^(J−1) f _(s1) h _(j−s)(θ)  (59)where, for simplicity, our notation for f_(j2) does not account for itsdependence on θ. Moreover, the corresponding least-squares objectivefunction, L_(θ), has the same form as L (see equation (4))

$\begin{matrix}{{{L_{\theta}\left( {c,f_{1},f_{3},\theta} \right)} = {\sum\limits_{i = 1}^{I}\;{\sum\limits_{j = 0}^{J - 1}\;\left( {a_{ij} - {\sum\limits_{k = 1}^{K}\;{c_{ik}f_{jk}}}} \right)^{2}}}},} & (60)\end{matrix}$except it depends on θ because f_(j2)=f_(j1)*h_(j)(θ). To insure thath_(j)(θ) is a nonnegative function, the following feasible set for θ ischosen:D _(θ)

{g,p,m,d: g≥0,0≤p≤1,m=0, 1, 2, . . . , d=0, 1, 2, . . . }.  (61)Hence, for θ∈D_(θ), ƒ_(j2) is a nonnegative function provided ƒ_(j1) isa nonnegative function.

Convolution Least-Squares Algorithm

Given the similarity of L_(θ) and L, it follows that the correspondingmajorizing functions for L_(θ) are

$\begin{matrix}{{q_{\theta}\left( {c,c^{(n)},f_{1},f_{3},\theta} \right)}\overset{\Delta}{=}{\sum\limits_{i = 1}^{I}\;{\sum\limits_{j = 0}^{J - 1}\;\left\{ {a_{ij}^{2} - {2\; a_{ij}{\sum\limits_{k = 1}^{K}\;{c_{ik}f_{jk}}}} + {{{\hat{a}}_{ij}\left( {c^{(n)},f} \right)}{\sum\limits_{k = 1}^{K}\;{\frac{c_{ik}^{2}}{c_{ik}^{(n)}}f_{jk}}}}} \right\}}}} & (62) \\{{r_{\theta}\left( {f_{1},f_{3},f_{1}^{(n)},f_{3}^{(n)},c,\theta} \right)}\overset{\Delta}{=}{\sum\limits_{i = 1}^{I}\;{\sum\limits_{j = 0}^{J - 1}\;\left\{ {a_{ij}^{2} - {2\; a_{ij}{\sum\limits_{k = 1}^{K}\;{c_{ik}f_{jk}}}} + {{{\hat{a}}_{ij}\left( {c,f^{(n)}} \right)}{\sum\limits_{k = 1}^{K}\;{\frac{f_{jk}^{2}}{f_{jk}^{(n)}}c_{ik}}}}} \right\}}}} & (63)\end{matrix}$where, again, fj2=fj1*hj(θ). By construction, for all c, c^((n))∈D_(c),f₁, f₃≥0, and θ∈D_(θ),q _(θ)(c,c ^((n)) ,f ₁ ,f ₃,θ)≥L _(θ)(c,f ₁ ,f ₃,θ)  (64)q _(θ)(c,c,f ₁ ,f ₃,θ)=L _(θ)(c,f ₁ ,f ₃,θ)  (65)r _(θ)(f ₁ ,f ₃ ,f ₁ ^((n)) ,f ₃ ^((n)) ,c,θ)≥L _(θ)(c,f ₁ ,f ₃,θ)  (66)r _(θ)(f ₁ ,f ₃ ,f ₁ ,f ₃ ,c,θ)=L _(θ)(c,f ₁ ,f ₃,θ)  (67)Using the MM methodology, the updates for estimating the factor weightsand tissue curves, given the current parameters estimates, are given by

$\begin{matrix}{c^{({n + 1})} = {\arg{\min\limits_{c \geq 0}{q_{\theta}\left( {c,c^{(n)},f_{1}^{(n)},f_{3}^{(n)},\theta^{(n)}} \right)}}}} & (68) \\{\left( {f_{1}^{({n + 1})},f_{3}^{({n + 1})}} \right) = {\arg{\min\limits_{f_{1},{f_{3} \geq 0}}{r_{\theta}\left( {f_{1},f_{3},f_{1}^{(n)},f_{3}^{(n)},c^{({n + 1})},\theta^{(n)}} \right)}}}} & (69) \\{\theta^{({n + 1})} = {\arg{\min\limits_{\theta \in D_{\theta}}{L_{\theta}\left( {c^{({n + 1})},f_{1}^{({n + 1})},f_{3}^{({n + 1})},\theta} \right)}}}} & (70)\end{matrix}$

The optimization problem in equation (69) is not straightforward becausethe objective function in (69) has “cross terms” of the formf_(j1)f_(j′1), j≠j′. These cross terms are due to the f_(j2) ² term inequation (63). Thus, it will be beneficial to construct a majorizingfunction for f_(j2) ². It is straightforward to show that the steps usedto derive the inequality in equation (22) can be repeated to get thefollowing inequality

$\begin{matrix}{{{f_{j\; 2}^{2} \leq {\upsilon\left( {f_{1},f_{1}^{(n)},\theta} \right)}}\overset{\Delta}{=}{\left\lbrack {\sum\limits_{s^{\prime} = 0}^{J - 1}\;{f_{s^{\prime}1}^{(n)}{h_{j - s^{\prime}}(\theta)}}} \right\rbrack{\sum\limits_{s = 0}^{J - 1}\;{{h_{j - s}(\theta)}\frac{f_{s\; 1}^{2}}{f_{s\; 1}}}}}},} & (71)\end{matrix}$which has the property that υ(f₁, f₁, θ)=f _(j2) ² for all f₁≥0 andθ∈D_(θ). Now, equation (71) is substituted into equation (63) to get amajorizing function that satisfies equations (66) and (67) and can beeasily minimized with respect to f₁

$\begin{matrix}{\mspace{79mu}{{\overset{\_}{r}}_{\theta}\left( {f_{1},f_{3},f_{1}^{(n)},f_{3}^{(n)},c,\theta} \right)}} & (72) \\{{\overset{\Delta}{=}{\sum\limits_{i = 1}^{I}\;{\sum\limits_{j = 0}^{J - 1}\;\left\{ {a_{ij}^{2} - {2\; a_{ij}{\sum\limits_{k = 1}^{K}\;{c_{ik}f_{jk}}}} + {{{\hat{a}}_{ij}\left( {c,f^{(n)}} \right)}\left( {{c_{i\; 1}\;\frac{f_{j\; 1}^{2}}{f_{j\; 1}^{(n)}}} + {c_{i\; 1}{\sum\limits_{s = 0}^{J - 1}\;{{h_{j - s}(\theta)}\frac{f_{s\; 1}^{2}}{f_{s\; 1}}}}} + {c_{i\; 3}\frac{f_{j\; 3}^{2}}{f_{j\; 3}^{(n)}}}} \right)}} \right\}}}},{where}} & (73) \\{\mspace{79mu}{f_{j\; 2}^{(n)} = {{f_{j\; 1}^{(n)}*{h_{j}\left( \theta^{(n)} \right)}} = {\sum_{s = 0}^{J - 1}{f_{s\; 1}^{(n)}{h_{j - s}\left( \theta^{(n)} \right)}}}}}} & (74)\end{matrix}$

Comparing equations (24) with (62), and (25) with (63), it is evidentthat the update for the factor weights c (see equation (68)) and thethird factor curve f₃ (see equation (69)) are given by equations (28)and (29), respectively, except that now f_(j2) ^((n)) is given byequation (74). In order to get the update for the right ventricle f₁(see equation (69)), we take the partial derivative of r _(θ)((f₁, f₃,f₁ ^((n)), f₃ ^((n)), c^((n+1)), θ^((n))) with respect to f₁ and set theresult to zero. From straightforward calculations, the derivative of r_(θ)((f₁, f₃, f₁ ^((n)), f₃ ^((n)) c^((n+1)), θ^((n))) with respect tof_(j′1) is

$\begin{matrix}{\frac{\partial\;{{\overset{\_}{r}}_{\theta}\left( {f_{1},f_{3},f_{1}^{(n)},f_{3}^{(n)},c^{({n + 1})},\theta^{(n)}} \right.}}{\partial f_{j^{\prime}1}} = {{{- 2}{\sum\limits_{i = 1}^{I}\;{\sum\limits_{j = 0}^{J - 1}{c_{i\; 1}a_{{ij}^{\prime}}}}}} - {2{\sum\limits_{i = 1}^{I}\;{\sum\limits_{j = 0}^{J - 1}\;{c_{i\; 2}a_{ij}{h_{j - j^{\prime}}\left( \theta^{(n)} \right)}}}}} + {2\frac{f_{j^{\prime}1}}{f_{j},_{1}^{(n)}}{\sum\limits_{i = 1}^{I}\;{c_{i\; 1}{{\hat{a}}_{{ij}^{\prime}}\left( {c^{({n + 1})},f^{(n)}} \right)}}}} + {2\frac{f_{j^{\prime}1}}{f_{j^{\prime}1}^{(n)}}{\sum\limits_{i = 1}^{I}{\sum\limits_{j = 0}^{J - 1}\;{c_{i\; 2}{{\hat{a}}_{{ij}^{\prime}}\left( {c^{({n + 1})},f^{(n)}} \right)}{h_{j - j^{\prime}}\left( \theta^{(n)} \right)}}}}}}} & (75) \\{= {{{- 2}{\sum\limits_{i = 1}^{I}\;{c_{i\; 1}{\hat{a}}_{ij}}}} - {2{\sum\limits_{i = 1}^{I}\;{c_{i\; 2}\left\lbrack {a_{ij}*{h_{- j^{\prime}}\left( \theta^{(n)} \right)}} \right\rbrack}}} + {2\frac{f_{j^{\prime}1}}{f_{j^{\prime}1}^{(n)}}{\sum\limits_{i = 1}^{I}\;{c_{i\; 1}{{\hat{a}}_{{ij}^{\prime}}\left( {c^{({n + 1})},f^{(n)}} \right)}}}} + {2\frac{f_{j^{\prime}1}}{f_{j^{\prime}1}^{(n)}}{\sum\limits_{i = 1}^{I}\;{{c_{i\; 2}\left\lbrack {{{\hat{a}}_{{ij}^{\prime}}\left( {c^{({n + 1})},f^{(n)}} \right)}*{h_{- j^{\prime}}\left( \theta^{(n)} \right)}} \right\rbrack}.}}}}} & (76)\end{matrix}$

Now, equation (76) is set to zero to get the desired update for theright ventricle tissue curve, which is

$\begin{matrix}\begin{matrix}{{f_{j\; 1}^{({n + 1})} = {f_{j\; 1}^{(n)}\frac{{\sum_{i = 1}^{I}\;{c_{i\; 1}a_{i\; j}}} + {\sum_{i = 1}^{I}\;{c_{i\; 2}\left\lbrack {a_{ij}*{h_{- j}\left( \theta^{(n)} \right)}} \right\rbrack}}}{{\sum_{i = 1}^{I}\;{c_{i\; 1}{{\hat{a}}_{ij}\begin{pmatrix}{c^{({n + 1})},} \\f^{(n)}\end{pmatrix}}}} + {\sum_{i = 1}^{I}\;{c_{i\; 2}\begin{bmatrix}{{{\hat{a}}_{ij}\begin{pmatrix}{c^{({n + 1})},} \\f^{(n)}\end{pmatrix}}*} \\{h_{- j}\left( \theta^{(n)} \right)}\end{bmatrix}}}}}},j} \\{{= 0},1,\ldots\mspace{11mu},{J - 1.}}\end{matrix} & (77)\end{matrix}$When knowledge about the maximum value of the right ventricle tissuecurve is incorporated, as well as the penalty function Λ, then theupdate for the right ventricle tissue curve is

$\begin{matrix}\begin{matrix}{{f_{j\; 1}^{({n + 1})} = \left\lbrack {f_{j\; 1}^{(n)}\frac{{\sum_{i = 1}^{I}\;{c_{i\; 1}a_{i\; j}}} + {\sum_{i = 1}^{I}\;{c_{i\; 2}\left\lbrack {a_{ij}*{h_{- j}\left( \theta^{(n)} \right)}} \right\rbrack}}}{{\sum_{i = 1}^{I}\;{c_{i\; 1}{{\hat{a}}_{ij}\begin{pmatrix}{c^{({n + 1})},} \\f^{(n)}\end{pmatrix}}}} + {\sum_{i}^{I}\;{c_{i\; 2}\begin{bmatrix}{{{\hat{a}}_{ij}\begin{pmatrix}{c^{({n + 1})},} \\f^{(n)}\end{pmatrix}}*} \\{h_{- j}\left( \theta^{(n)} \right)}\end{bmatrix}}}}} \right\rbrack_{\mu_{k} - \epsilon}^{\mu_{k} + \epsilon}},j} \\{{= 0},1,\ldots\mspace{11mu},{J - 1.}}\end{matrix} & (78)\end{matrix}$

With the update for the right ventricle curve in hand (i.e., update inequations (77) or (78)), the algorithm for minimizing the least-squaresfunction L_(θ), is proposed, which is referred to herein as theconvolution least-squares (CLA) algorithm. First, let d_(min) andd_(max) denote the minimum and maximum delay considered, respectively,and m_(max) denote the maximum value considered for the parameter m. TheCLA algorithm follows:

get initial estimates: f₁ ⁽⁰⁾, f₃ ⁽⁰⁾, g⁽⁰⁾, and p⁽⁰⁾ for {circumflexover (d)} = d_(min), d_(min) + 1, ... , d_(max),  for {circumflex over(m)} = 1,2, ... , m_(max),   for n=0,1,2, . . .    let θ^((n)) =[g^((n)) p^((n)) {circumflex over (m)} {circumflex over (d)}]    getc^((n+1)) using (28) with f_(j2) ^((n)) given by (74)    get f₁ ^((n+1))using (77)    get f₃ ^((n+1)) using (29) with f_(j2) ^((n)) given by(74) (79)    $\left( {g^{({n + 1})},p^{({n + 1})}} \right) = {\arg{\min\limits_{{\theta \geq 0},{0 \leq p \leq 1}}{L_{\theta}\left( {c^{({n + 1})},f_{1}^{({n + 1})},f_{3}^{({n + 1})},g,p,\hat{m},\hat{d}} \right)}}}$  end   store current parameter estimates and corresponding value forL_(θ)  end end

The desired parameter estimates are the estimates from the CLA algorithmthat produce the smallest least-squares error. It should be noted thatthe objective function L_(θ) decreases monotonically with increasingiterations.

The problem in equation (79) can be solved using the following iterativecoordinate descent method. It will be convenient to express h_(j)(θ) ash_(j)(θ)=gh _(j)(p, m, d), where h _(j)(p, m, d) is the inverseZ-transform of

$\begin{matrix}{{\overset{\_}{H}(z)} = {\left( \frac{1 + z^{- 1}}{1 - {pz}^{- 1}} \right)^{m}{z^{- d}.}}} & (80)\end{matrix}$With this notation, the objective function in (79) can be written as

$\begin{matrix}{{{L_{\theta}\left( {c^{({n + 1})},f_{1}^{({n + 1})},f_{3}^{({n + 1})},g,p,m,\hat{d}} \right)} = {{\alpha^{({n + 1})}g^{2}} - {2\;\beta^{({n + 1})}g} + \gamma^{({n + 1})}}},\mspace{79mu}{where}} & (81) \\{{\alpha^{({n + 1})}\left( {c^{({n + 1})},f_{1}^{({n + 1})},f_{3}^{({n + 1})}} \right)}\overset{\Delta}{=}{\sum\limits_{i = 1}^{I}\;{\sum\limits_{j = 0}^{J - 1}\;{\left\lbrack {a_{ij} - \left( {{c_{i\; 1}^{({n + 1})}f_{j\; 1}^{({n + 1})}} + {c_{i\; 3}^{({n + 1})}f_{j\; 3}^{({n + 1})}}} \right)} \right\rbrack{\beta^{({n + 1})}\left( {c^{({n + 1})},f_{1}^{({n + 1})},f_{3}^{({n + 1})},p,\hat{m},\hat{d}} \right)}}}}} & (82) \\{\overset{\Delta}{=}{\sum\limits_{i = 1}^{I}\;{\sum\limits_{j = 0}^{J - 1}\;{\left\lbrack {a_{ij} - \left( {{c_{i\; 1}^{({n + 1})}f_{j\; 1}^{({n + 1})}} + {c_{i\; 3}^{({n + 1})}f_{j\; 3}^{({n + 1})}}} \right)} \right\rbrack \times {c_{i\; 2}^{({n + 1})}\left( {f_{j\; 1}^{({n + 1})}*{{\hat{h}}_{j}\left( {\hat{p},\hat{m},\hat{d}} \right)}} \right)}}}}} & (83) \\{{\gamma^{({n + 1})}\left( {c^{({n + 1})},f_{1}^{({n + 1})},f_{3}^{({n + 1})},p,\hat{m},\hat{d}} \right)}\overset{\Delta}{=}{\sum\limits_{i = 1}^{I}\;{\sum\limits_{j = 0}^{J - 1}\;{\left\lbrack {c_{i\; 2}^{({n + 1})}\left( {f_{j\; 1}^{({n + 1})}*{{\hat{h}}_{j}\left( {\hat{p},\hat{m},\hat{d}} \right)}} \right)} \right\rbrack^{2}.}}}} & (84)\end{matrix}$Thus, the unconstrained minimizer of L_(θ)(c^((n+1)), f₁ ^((n+1)), f₃^((n+1)), g, {circumflex over (p)}, {circumflex over (m)}, {circumflexover (d)}) with respect to g is simply

$\begin{matrix}{{\hat{g} = \frac{\beta^{({n + 1})}\left( {c^{({n + 1})},f_{1}^{({n + 1})},f_{3}^{({n + 1})},\hat{p},\hat{m},\hat{d}} \right)}{\alpha^{({n + 1})}\left( {c^{({n + 1})},f_{1}^{({n + 1})},f_{3}^{({n + 1})},\hat{p},\hat{m},\hat{d}} \right)}},} & (84)\end{matrix}$where {circumflex over (p)} denotes the current estimate for the filterparameter p. With the above result, the steps for solving theoptimization problem of equation (79) is given using the coordinatedescent method: Let p^(old)=p^((n)) and [x]+=max(0, x)

-   IC Step 1:

$\begin{matrix}{g^{new} = {\arg{\min\limits_{g \geq 0}{L_{\theta}\left( {c^{({n + 1})},f_{1}^{({n + 1})},f_{3}^{({n + 1})},p^{old},\hat{m},\hat{d}} \right)}}}} & {{~~~~~~~~~~~~~~~~~~~~~~~~}(84)} \\{{= \left\lbrack \frac{\beta^{({n + 1})}\left( {c^{({n + 1})},f_{1}^{({n + 1})},f_{3}^{({n + 1})},p^{old},\hat{m},\hat{d}} \right)}{\alpha^{({n + 1})}\left( {c^{({n + 1})},f_{1}^{({n + 1})},f_{3}^{({n + 1})},p^{old},\hat{m},\hat{d}} \right)} \right\rbrack_{+}},} & {(87)}\end{matrix}$

-   IC Step 2: Using a 1D line search, such as the golden section method    [17], get p^(new):

$\begin{matrix}{p^{new} = {\arg{\min\limits_{0 \leq p \leq 1}{L_{\theta}\left( {c^{({n + 1})},f_{1}^{({n + 1})},f_{3}^{({n + 1})},g^{new},p,\hat{m},\hat{d}} \right)}}}} & (88)\end{matrix}$

-   IC Step 3: Set p^(old)=p^(new) and repeat the above two steps, for a    set number of iterations or until some desired stopping criterion is    met.-   IC Step 4: Set g^((n+1))=g^(new) and p^((n+1))=p^(new).

So configured, the CLA algorithm is monotonic and guaranteed to producenonnegative estimates.

Initialization

User-dependent and user-independent methods are available fordetermining initial estimates for the right and left ventricle tissuescurves, and myocardial tissue curve. The initial factor weights aretypically chosen to be uniform in the sense that c_(ik)=⅓ for all i, k.Regarding the range of values for the delay, suitable values for theminimum and maximum delay between the right and left ventricle tissuecurves are available from the literature in human physiology.

It is expected that a suitable maximum value for the filter order,m_(max), can be determined through experiments. Consequently, thefollowing addresses computing initial estimates for the filterparameters g and p.

Because F₂(z)=H(z)F₁(z) from equation (59), it follows that F₂⁽⁰⁾(z)≈H(z)F₁ ⁽⁰⁾(z), where f_(1j) ⁽⁰⁾ and f_(2j) ⁽⁰⁾ denote the initialright and left ventricle tissue curves, and H(z) is given by equation(58). Therefore,

$\begin{matrix}{{F_{2}^{(0)}(z)} \approx {{H(z)}{F_{1}^{(0)}(z)}}} & (89) \\{{F_{2}^{(0)}(z)} \approx {{g\left( \frac{1 + z^{- 1}}{1 - {pz}^{- 1}} \right)}^{m}z^{- d}{F_{1}^{(0)}(z)}}} & (90) \\{{{F_{2}^{(0)}(z)} \approx {g\frac{1}{\left( {1 - {pz}^{- 1}} \right)^{m}}z^{- d}{{\overset{\_}{F}}_{1}^{(0)}(z)}}}{where}} & (91) \\{{{\overset{\_}{F}}_{1}^{(0)}(z)} = {{F_{1}^{(0)}(z)}\left( {1 + z^{- 1}} \right)^{m}{z^{- d}.}}} & (92)\end{matrix}$

Thus, the initial left ventricle tissue curve is assumed to beapproximately equal to the response of a certain mth-order, all-polefilter to the initial right ventricle tissue curve. This observation andequation (91) forms the basis of the following method for obtaininginitial estimates for g and p:

 for {circumflex over (d)} = d_(min), d_(min) + 1, ... , d_(max),   form = 1, 2, ... , m_(max), Initial g, p - Step 1: Using an input-outputsystem identification method, such as the Steiglitz-McBride algorithm,find the parameters b₀, a₁, a₂, ... , a_(m) that provide the bestleast-squares fit for the following model: (93)${F_{2}^{(0)}(z)} = {b_{0}\frac{1}{\left( {1 + {a_{1}z^{- 1}} + {a_{2}z^{- 2}} + \ldots + {a_{m}z^{- m}}} \right)}{{\overset{\_}{F}}_{1}^{(0)}(z)}}$Initial g, p - Step 2: Let {circumflex over (b)}₀, â₁, â₂, . . . , â_(m)be the resulting estimates for b₀, a₁, a₂, ... , a_(m), respectively,from the first step. Comparing equations (91) and (93), an estimate forg is ĝ = {circumflex over (b)}₀. Moreover, an estimate for p can beobtained by determining the best least-squares fit between (1 -pz⁻¹)^(m) and (1 + a₁z⁻¹ +a₂z⁻² + ... + a_(m)z^(-m)). For Example, whenm = 2, an estimate for p is: (94)${\hat{p} = {{\arg{\min\limits_{0 \leq p \leq 1}\left( {{{- 2}p} - {\hat{a}}_{1}} \right)^{2}}} + \left( {p^{2} - {\hat{a}}_{2}} \right)^{2}}},$Initial g, p - Step 3: Store the current parameter estimates {ĝ,{circumflex over (p)}, {circumflex over (m)}, {circumflex over (d)}} andthe corresponding least-squares error (95)${{L_{init}\left( {\hat{g},\hat{p},\hat{m},\hat{d}} \right)}\overset{\Delta}{=}{\sum\limits_{j = 0}^{J - 1}\left( {f_{j\; 2}^{(0)} - {{{\overset{\_}{f}}_{j\; 1}^{(0)}\left( {\hat{m},\hat{d}} \right)}*{{\overset{\sim}{h}}_{j}\left( {\hat{g},\hat{p},\hat{m}} \right)}}} \right)^{2}}},$where f _(j1) ⁽⁰⁾ ({circumflex over (m)}, {circumflex over (d)}) and{tilde over (h)}_(j)(g, p, m) are, respectively, the inverse        Z-transforms of F ₁ ⁽⁰⁾ (z) (see (92)) and (96)${\overset{\sim}{H}(z)}\overset{\Delta}{=}{g{\frac{1}{\left( {1 - {pz}^{- 1}} \right)^{m}}.}}$  end  end

The set of parameters that produce the smallest least-squares error arethe desired initial estimates for the parameters g, p, m, and d, whichwe denote respectively as g⁽⁰⁾, p⁽⁰⁾, m⁽⁰⁾, and d⁽⁰⁾.

The initial estimates for the parameters m and d are not explicitly usedin the CLA algorithm. However, they could be used to reduce the range ofvalues considered for the parameters m and d. For example, the “forloop” for the delay in the CLA algorithm could be instead: for{circumflex over (d)}=d⁽⁰⁾−Δ: d⁽⁰⁾+Δ, where Δ>0, would be an integerchosen by the user.

III. An Approach without Assumption of a Mathematical RelationshipBetween Right Ventricle and Left Ventricle Curves

In a further approach to estimating MBF, we first incorporate theHutchins pharmacological kinetic model into the standard FADS model in away to not assume that the RV and LV tissue curves obey any particularmathematical relationship unlike the methods described above. Theimproved model, which we refer to as the second pharmacological kineticsbased FADS (K-FADS-II) model, provides another way for estimating voxeldependent myocardium tissue curves that are physiologically meaningful.In our next step, we perform a discretization that transforms thecontinuous-time K-FADS-II model into a discrete-time K-FADS-II model. Itshould be noted that there is a simple relationship between thediscrete-time and continuous-time K-FADS-II parameters. Lastly, wedevelop an algorithm, which we call the Improved Voxel-ResolutionMyocardial Blood Flow (IV-MBF) algorithm that estimates the parametersof the discrete-time K-FADS-II model by iteratively minimizing a certainleast-squares (LS) objective function. The desired MBF estimates arecomputed in a straightforward manner from the estimated discrete-timeK-FADS-II model parameters. The IV-MBF algorithm was evaluatedsubjectively and objectively through experiments conducted usingsynthetic data and patient data, and found to be accurate and stable inthe test scenarios we considered. Hence, we believe the proposed methodis feasible for determining physiologically meaningful estimates of MBF.

To avoid confusion in view of the discussion above regarding otherapproaches, we start again from the fundamentals and restart thenumbering of equations. The kinetic behavior of ammonia in themyocardium is modeled by the compartment model shown in FIG. 1 and thefollowing differential equations (note: we assume k₄=0).

$\begin{matrix}{\frac{d\;{C_{F}(t)}}{d\; t} = {{k_{1}{C_{P}(t)}} - {\left( {k_{2} + k_{3}} \right){C_{F}(t)}}}} & (1) \\{{\frac{d\;{C_{T}(t)}}{d\; t} = {k_{3}{C_{F}(t)}}},} & (2)\end{matrix}$where C_(P) is the ammonia concentration in blood plasma, C_(F) is theconcentration of the free (i.e., unbound) ammonia, and C_(T) is theconcentration of the trapped (i.e., bound) ammonia.

In the standard FADS model, it is assumed that a_(ij), the activity invoxel i of frame j, is a linear combination of K primary factor curves

$\begin{matrix}{a_{i,j} = {\sum\limits_{k = 1}^{K}\;{c_{ik}f_{kj}}}} & (3)\end{matrix}$where {f_(kj)} are the values of the factor curves and the coefficients{c_(ik)} are the factor weights. The primary factor curves forconventional MBF estimation applications are the right ventricle (RV),left ventricle (LV), and myocardial tissue curves, which model theammonia concentration as a function of time in the RV, LV, andmyocardium, respectively. The mathematical task is to find both thefactor curves and weights so that the linear combination of factorcurves for every voxel in the myocardium matches the correspondingmeasured time activity curve (TAC) as close as possible. This problem isconstrained by requiring that the factor curves and weights all benonnegative.

Other approaches include spillover correction because of the finiteresolution of the scanner and because the myocardium is moving duringthe scan. In one known approach, factor analysis was used to obtain, intheory, spillover independent time activity curves of the RV, LV, andmyocardial blood tissue. By using curves generated from factor analysis,the spillover component in the model can be eliminated in theory.However, it is noted that factor analysis does not correct for the undermeasurement due to the partial volume effect. Correcting errors due tothe partial volume effect would require the use of a contrast recoverycoefficient.

A. K-Fads-II Model

In this section, we first present a model that combines Hutchins'pharmacological kinetic model with the standard factor analysis model(see equation (3)). The resulting model is an improvement over anearlier version discussed above, so we call it the secondpharmacological kinetics based FADS model (K-FADS-II). Next, we discussa discretization that transforms the continuous-time K-FADS-II modelinto a discrete-time K-FADS-II model. Specifically, the discrete-timeK-FADS-II model is obtained by applying the bilinear transform to thecontinuous-time K-FADS-II model. It should be noted that there is asimple relationship between the discrete-time and continuous-timeK-FADS-II parameters.

1. Continuous-Time Model for Activity in the Myocardium

The unknown ammonia concentrations in the RV and LV (i.e., C_(P)) areassumed to be spatially independent, whereas the unknown concentrationsof the free ammonia (i.e., C_(F)) and trapped ammonia (i.e., C_(T)) areassumed to spatially dependent. With these assumptions in mind, we letthe RV and LV tissue curves be denoted by the continuous-time functionsƒ(t) and g(t), respectively, and the free and trapped ammoniaconcentrations in the i^(th) voxel within the myocardium (i.e., thei^(th) myocardium voxel) be denoted by g_(i,F)(t), and g_(i,T)(t),respectively.

Under the above assumptions it follows that the rate constants inHutchin's model are spatially dependent. After applying the Laplacetransform to voxel dependent versions of equations (1) and (2), we getthe following expressions for the concentrations of the free and trappedammonia in the i^(th) voxel within the myocardiumg _(i,F)(t)=k _(i,1) g(t)*exp(− k _(i) t)u(t)  (4)g _(i,T)(t)=k _(i,3) g _(i,F)(t)*u(t)  (5)where k _(i)

k_(i,2)+k_(i,3), I is the number of voxels within the myocardium, and *denotes convolution operator. We note that k_(1,1), k_(2,1), . . . ,k_(I,1) are the desired MBF parameters. In keeping with the assumptionsbehind (3), we assume the true TAC for the i^(th) voxel in themyocardium is a superposition of the ammonia activity in the RV and LV,and free and trapped ammonia activity in the myocardial tissuea _(i)(t)=c _(i,1) f(t)+c _(i,2) g(t)+c _(i,3) g _(i,F)(t)+c _(i,4) g_(i,T)(t),i=1,2 . . . ,I.  (6)The first term in (6) is identified as the amount of spillover from theRV, and the second term is a combination of the ammonia activity inblood vessels within the myocardium and spillover from the LV. Morespecifically, the constant c_(i,1) accounts for the amount of theammonia activity in voxel i that is due to the blood plasma in the RV.Further, the constant c_(i,2) accounts for the amount of the ammoniaactivity in voxel i that is due to the blood plasma in the LV (i.e., LVspill over) and blood plasma in the blood vessels of the myocardium. Thethird and fourth terms in (6) model the free and trapped ammoniaactivity in the myocardium, respectively. The coefficients c_(i,3) andc_(i,4) represent the fractional volume of voxel i that can be occupiedby the ammonia activity in either the free or trapped states,respectively. Given that the free space for water in myocardial tissueis approximately 80%, we assume that c_(i,3)=c_(i,4)=0.8.

2. Discrete-Time Model for Activity in Myocardium

When a dynamic PET scan of the heart is taken, the images are separatedinto three sets of voxels that lie in the RV, LV, and myocardium such asthat described by V. Appia, B. Ganapathy, A. Yezzi, and T. Faber in“Localized principal component analysis based curve evolution: A divideand conquer approach,” 2011 International Conference on Computer Vision,pp. 1981-1986, 2011, which is incorporated by reference. Let themeasured TACs associated with the set of voxels in the RV, LV, andmyocardium be referred to as the RV TACs, LV TACs, and myocardium TACs,respectively. To model the underlying sampling process of dynamic PETscanning protocols, we let d_(i)[n], r[n], and l[n] represent thediscrete-time signals obtained by sampling a_(i)(t), f(t), and g(t),respectivelyd _(i) [n]

a _(i)(nT _(s)), n=0, 1, . . . , N−1, i=1, 2, . . . , I  (7)r[n]

f(nT _(s)),n=0,1, . . . ,N−1,  (8)l[n]

g(nT _(s)),n=0,1, . . . ,N−1,  (9)where

$f_{s}\overset{\Delta}{=}\frac{1}{T}$is the sampling rate and T is the sampling interval. We will refer tothe quantities d_(i)[n] as the i^(th) sampled myocardium TAC, and thediscrete time signals r[n] and l[n], which are unknown, as the sampledRV tissue curve and sampled LV tissue curve, respectively. Inapplications where the scan durations for dynamic sequence PET protocolsare not uniform, the assumption that the TACs are sampled uniformly isinappropriate. However, uniform samples of the TACs can be obtained fromnon-uniform samples of the TACs via a suitable interpolation using knownmethods.

It is of interest to determine a discrete-time system with the propertythat its response to the sampled RV and LV tissue curves closelyapproximates the i^(th) sampled myocardium TAC. The bilineartransformation is a popular way to transform a linear time-invariantcontinuous-time system into a linear time-invariant discrete-timesystem. Specifically, the system function of a continuous-time system,H(s), can be converted into the system function of a discrete-timesystem H_(d)(z) using

$\begin{matrix}{{H_{d}(z)} = \left. {H(s)} \right|_{s = {\frac{2}{T_{s}}\frac{1 - z^{- 1}}{1 + z^{- 1}}}}} & (10)\end{matrix}$In our analysis, we make use of the z-transform. Let x[n] be anarbitrary discrete-time sequence. The z-transform of x[n] is defined tobe X(z)=Σ_(n=−∞) ^(∞)x[n]z^(−n).

Taking the Laplace transform of (6), we get the following relationship

$\begin{matrix}{{{A_{i}(s)} = {{c_{i,1}{F(s)}} + {c_{i,2}{G(s)}} + {{G(s)}{H_{i,F}(s)}} + {{G(s)}{H_{i,T}(s)}}}}{where}} & (11) \\{{H_{i,F}(s)}\overset{\Delta}{=}{c_{i,3}\frac{k_{i,1}}{s + k_{i}}}} & (12) \\{{H_{i,T}(s)}\overset{\Delta}{=}{c_{i,4}\frac{k_{i,1}}{s + k_{i}}{\frac{k_{i,3}}{s}.}}} & (13)\end{matrix}$After sampling a_(i)(t), and applying the bilinear transformation to thesystem functions H_(i,F)(s) and H_(i,T)(s), we get the desireddiscrete-time model

$\begin{matrix}{{{D_{i}(z)} = {{c_{i,1}{R(z)}} + {c_{i,2}{L(z)}} + {c_{i,3}{k_{i,1}\left( {{2\text{/}T_{s}} + k_{i}} \right)}^{- 1}{L(z)}\frac{1 + z^{- 1}}{1 - {p_{i}z^{- 1}}}} + {c_{i,4}k_{i,1}{k_{i,3}\left( {{2\text{/}T_{s}} + k_{i}} \right)}^{- 1}\left( {2\text{/}T_{s}} \right)^{- 1}{L(z)}\frac{1 + z^{- 1}}{1 - {p_{i}z^{- 1}}}\frac{1 + z^{- 1}}{1 - z^{- 1}}}}}\mspace{79mu}{{{where}\mspace{14mu} p_{i}}\overset{\Delta}{=}{\left( {{2\text{/}T_{s}} - k_{i}} \right){\left( {{2\text{/}T_{s}} + k_{i}} \right)^{- 1}.}}}} & (14)\end{matrix}$

The bilinear transformation has the property that it maps a stablecontinuous-time system into a stable discrete-time system. Moreover, thebilinear transformation avoids the problem of aliasing by mapping the jΩaxis into the unit circle of the complex plane. However, frequencywarping occurs as a result of mapping the entire jΩ axis into the unitcircle. Note, the frequency warping problem can be ameliorated bychoosing a sufficiently high sampling rate f_(s).

From equation (14), it follows that the i^(th) myocardium TAC, which isnoisy in practice, can be expressed asd _(i)(n)=c _(i,1) r(n)+c _(i,2) l(n)+k _(i,1)(2/T _(s) +k _(i))⁻¹l[n]*h _(F) [n;p _(i) ]+k _(i,1) k _(i,3)(2/T _(s) +k _(i))⁻¹(2/T_(s))⁻¹ l[n]*h _(T) [n;p _(i) ]+e[n]  (15)where e[n] is the noise and * denotes discrete-time convolution. Also,

$\begin{matrix}{{h_{F}\left\lbrack {n;\alpha} \right\rbrack}\overset{\Delta}{=}{c_{i,3}{\mathcal{Z}^{- 1}\left\lbrack \frac{1 + z^{- 1}}{1 - {\alpha\; z^{- 1}}} \right\rbrack}}} & (16) \\{{h_{T}\left\lbrack {n;\alpha} \right\rbrack}\overset{\Delta}{=}{c_{i,4}{\mathcal{Z}^{- 1}\left\lbrack \frac{\left( {1 + z^{- 1}} \right)^{- 2}}{\left( {1 - {\alpha\; z^{- 1}}} \right)\left( {1 - z^{- 1}} \right)} \right\rbrack}}} & (17)\end{matrix}$where Z⁻¹ denotes the inverse z-transform. Given the data {d_(i)[n]},the problem is estimate the K-FADS-II model parameters:r

[r[0], r[1], . . . ,r[N−1]]  (18)l

[l[0], l[1], . . . ,l[N−1]]  (19)c

{c_(1,1), c_(2,1), . . . , c_(I,1), c_(1,2), c_(2,2), . . . ,c_(I,2)}  (20)k₁

[k_(1,1), k_(2,1), . . . , k_(I,1)]  (21)k₂

[k_(1,2), k_(2,2), . . . , k_(I,2)]  (22)k₃

[k_(1,3), k_(2,3), . . . , k_(I,3)]  (23)Recall the parameters k₁, are the desired MBF parameters. Thus, theparameters r, l, c, k₂, and k₃ can be viewed as nuisance parametersbecause they must be accounted for in the analysis even though they arenot of direct interest.

In the next section, we develop an algorithm for estimating the sampledRV and LV tissue curves and K-FADS-II parameters that we call theImproved Voxel Resolution MFB (IV-MBF) algorithm. The IV-MBF algorithmis based on a model that accounts for the fact that the shape of TACsdue to ischemic and normal tissue are different. In fact, the modelallows for tissue curves that represent free and trapped ammonia to bevoxel dependent and physiologically appropriate. By contrast, in thestandard FADS model, it is assumed that TACs in ischemic and normaltissue can be modeled as a linear combination of the same three factors.We believe the use of a more appropriate model offers the possibility ofmore accurate MBF estimates than available methods.

The IV-MBF algorithm performed well in a limited study using realpatient data. The results of the study suggest that the IV-MBF algorithmis robust and would perform well in practice, where MBF values due toischemic and normal tissue can vary over a wide range.

B. IV-MBF Algorithm

Consider the problem of minimizing a real valued function f with domainD∈R^(n). The majorize-minimize (MM) technique can be used to develop analgorithm that produces a sequence of iterates {x^((m))} such thatf(x^((m))) is a monotonically decreasing function such as thosedescribed by D. Hunter and K. Lange, “A Tutorial on MM Algorithms,” TheAmerican Statistician, vol. 58, pp. 30-37, 2004, which is incorporatedby reference. In this section, we first briefly introduce the MMtechnique and then we develop the IV-MBF algorithm by applying thistechnique to a certain LS objective function.

1. Review of the Majorize-Minimize Optimization Technique

A real valued function g with domain {(x₁, x₂): x₁, x₂∈D} is said tomajorize the function ƒ if the following conditions hold for all x, y∈D:g(x, y)≥ƒ(x)  (C1)g(x, x)=ƒ(x).  (C2)MM algorithms are a viable approach provided a majorizing function canbe found that is easier to minimize than the original objectivefunction. Assuming that g is a suitable majorizing function for ƒ, thecorresponding MM algorithm is

$\begin{matrix}{{x^{({k + 1})} = {\arg\mspace{14mu}{\min\limits_{x \in D}{g\left( {x,x^{(k)}} \right)}}}},} & (24)\end{matrix}$where x^((k)) is the current estimate for the minimizer of ƒ. It followsfrom (C1) and (C2) that the MM algorithm defined by equation (24)satisfies the monotonicity property [20]f(x ^((k+1)))≤g(x ^((k+1)) ,x ^((k)))≤g(x ^((k)) ,x ^((k)))=f(x^((k)))  (25)

2. IV-MBF Algorithm: LS Estimation of K-FADS-II Model Parameters

In our discussion we assume that suitable initial estimates areavailable. Later, in Section VII, we present a procedure for determininginitial estimates for the K-FADS-II model parameters.

1) Construction of a LS Objective Function L:

Because no statistical information about the noise is available, wepropose to estimate the K-FADS-II parameters using the LS estimationmethod. From equation (15), the LS objective function is given by

$\begin{matrix}{{L\left( {r,l,c,k_{1},k_{2},k_{3}} \right)}\overset{\Delta}{=}{\sum\limits_{i = 1}^{I}\;{\sum\limits_{n = 0}^{N - 1}\;\left\{ {{d_{i}\lbrack n\rbrack} - \left( {{c_{i,1}{r\lbrack n\rbrack}} + {c_{i,2}{l\lbrack n\rbrack}} + {{k_{i,1}\left( {{2\text{/}T_{s}} + k_{i}} \right)}^{- 1}{l\lbrack n\rbrack}*{h_{F}\left\lbrack {n;p_{i}} \right\rbrack}} + {k_{i,1}{k_{i,3}\left( {2\text{/}T_{s}} \right)}^{- 1}\left( {{2\text{/}T_{s}} + k_{1}} \right)^{- 1}{l\lbrack n\rbrack}*{h_{T}\left\lbrack {n;p_{i}} \right\rbrack}}} \right)} \right\}^{2}}}} & (26) \\{\mspace{205mu}{= {\sum\limits_{i = 0}^{I}\;{\sum\limits_{n = 0}^{N - 1}\;\left\{ {{d_{i}\lbrack n\rbrack} - \left( {{c_{i,1}{r\lbrack n\rbrack}} + {{l\lbrack n\rbrack}*{h\left\lbrack {n;\theta_{i,h}} \right\rbrack}}} \right)} \right\}^{2}}}}} & (27)\end{matrix}$where h[n; θ_(i,h)]

c_(i,2)δ[n]+k_(i,1)(2/T_(S)+k_(i))⁻¹h_(F)[n;p_(i)]+k_(i,1)k_(i,3)(2/T_(S))⁻¹(2/T_(S)+k_(i))⁻¹h_(T)[n; p_(i)] andθ_(i,h)

[c_(i,2), k_(i,1), k_(i,2), k_(i,3)]. It follows that a LS estimate ofthe K-FADS-II parameters is a solution to the following constrainedoptimization problem({circumflex over (r)},{circumflex over (l)},ĉ,{circumflex over (k)} ₁,{circumflex over (k)} ₂ ,{circumflex over (k)} ₃)=arg min L(r,l,c,k ₁,k ₂ ,k ₃) subject to r,l,c,k ₁ ,k ₂ ,k ₃≥0  (28)

b) Minimization of Least-Square Objective Function L:

We propose to minimize the LS objective function L by using the blockcoordinate descent method. In this method, the coordinates arepartitioned into a fixed number of blocks and, at each iteration, theobjective function is minimized with respect to one of the coordinateblocks while the remaining coordinates are fixed to their currentestimates. Let r^((m)), l^((m)), c^((m)), k₁ ^((m)), k₂ ^((m)), and k₃^((m)) denote the current estimates for r, l, c, k₁, k₂, and k₃,respectively. An outline of an algorithm that, in principle, could beused for minimizing L follows

Step 1 (29)$\left( {r^{({m + 1})},l^{({m + 1})}} \right)\overset{\Delta}{=}{\arg{\min\limits_{{r \geq 0},{l \geq 0}}{L\left( {r,l,c^{(m)},k_{1}^{(m)},k_{2}^{(m)},k_{3}^{(m)}} \right)}}}$Step 2 (30)$k_{2}^{({m + 1})}\overset{\Delta}{=}{\arg\mspace{11mu}{\min\limits_{k_{2} \geq 0}\;{L\left( {r^{({m + 1})},l^{({m + 1})},c^{(m)},k_{1}^{(m)},k_{2},k_{3}^{(m)}} \right)}}}$(31)$k_{3}^{({m + 1})}\overset{\Delta}{=}{\arg\mspace{11mu}{\min\limits_{k_{3} \geq 0}\;{L\left( {r^{({m + 1})},l^{({m + 1})},c^{(m)},k_{1}^{(m)},k_{2}^{({m + 1})},k_{3}} \right)}}}$Step 3 (32)$\left( {c^{({m + 1})},k_{1}^{({m + 1})}} \right)\overset{\Delta}{=}{\arg{\min\limits_{{c \geq 0},{k_{1} \geq 0}}{L\left( {r^{({m + 1})},l^{({m + 1})},c,k_{1},k_{2}^{({m + 1})},k_{3}^{({m + 1})}} \right)}}}$A simpler alternative to Step 1 is to use the MM technique describedabove to obtain iterates r^((m+1)) and l^((m+1)) that decrease theleast-squares objective L in the sense that, for all m=0, 1, 2, . . . ,L(r ^((m−1)) ,l ^((m+1)) ,c ^((m)) ,k ₁ ^((m)) ,k ₂ ^((m)) ,k ₃^((m)))≤L(r ^((m)) ,l ^((m)) ,c ^((m)) ,k ₁ ^((m)) ,k ₂ ^((m)) ,k ₃^((m)))  (33)Applying the same reasoning, we determine iterates c^((m+1)) and k₁^((m+1)) such thatL(r ^((m+1)) ,l ^((m+1)) ,c ^((m)) ,k ₁ ^((m+1)) ,k ₂ ^((m+1)) ,k ₃^((m+1)))≤L(r ^((m+1)) ,l ^((m+1)) ,c ^((m)) ,k ₁ ^((m)) ,k ₂ ^((m+1)),k ₃ ^((m+1)))  (34)Taking this approach, steps 1 and 3 of the proposed algorithm are now

[Step 1a] Find iterates r^((m+1)) and l^((m+1)) such thatL(r ^((m+1)) ,l ^((m+1)) ,c ^((m)) ,k ₁ ^((m)) ,k ₂ ^((m)) ,k ₃^((m)))≤L(r ^((m)) ,l ^((m)) ,c ^((m)) ,k ₁ ^((m)) ,k ₂ ^((m)) ,k ₃^((m)))  (35)

[Step 3a] Find iterates c^((m+1)) and k₁ ^((m+1)) such thatL(r ^((m+1)) ,l ^((m+1)) ,c ^((m+1)) ,k ₁ ^((m+1)) ,k ₂ ^((m+1)) ,k ₃^((m+1)))≤L(r ^((m+1)) ,l ^((m+1)) ,c ^((m)) ,k ₁ ^((m)) ,k ₂ ^((m+1)),k ₃ ^((m+1)))  (36)

Solution to Step 1a:

To construct a majorizing function for L at (r^(old), l^(old)), wherer^(old) and l^(old) are the current estimates for the sampled RV and LVtissue curves, respectively, we first construct a majorizing functionfor (c_(i,1)r[n]+l[n]*h[n; θ_(i,h)])² by exploiting the followingresult, which is proven in the Proof of Result 1 section below:

Result 1:

Let z₁[n], z₂[n], and w[n] be causal positive sequences (a sequence w[n]is said to be causal if w[n]=0 for n<0), where z₁[n] and z₂[n] are oflength N. A majorizing function for (z₁[n]+z₂[n]*w[n])² is

$\begin{matrix}{\left( {{z_{i}^{old}\lbrack n\rbrack} + {{z_{2}^{old}\lbrack n\rbrack}*{w\lbrack n\rbrack}}} \right)\left( {\frac{\left( {z_{1}\lbrack n\rbrack} \right)^{2}}{z_{1}^{old}\lbrack n\rbrack} + {\sum\limits_{k = 0}^{N - 1}\;{{w\left\lbrack {n - k} \right\rbrack}\frac{\left( {z_{2}\lbrack k\rbrack} \right)^{2}}{z_{2}^{old}\lbrack k\rbrack}}}} \right)} & (37)\end{matrix}$where z₁ ^(old)[n] and z₂ ^(old)[n] are causal positive sequences oflength N. From Result 1, it follows that a majorizing function for L at(r^(old), l^(old)) is

$\begin{matrix}{{q_{L}\left( {r,l,c,k_{1},k_{2},{k_{3};r^{old}},l^{old}} \right)} = {\sum\limits_{i = 1}^{I}\;{\sum\limits_{n = 0}^{N - 1}\;\left\{ {{d_{i}\lbrack n\rbrack}^{2} - {2{d_{i}\lbrack n\rbrack}\left( {{c_{i,1}{r\lbrack n\rbrack}} + {{l\lbrack n\rbrack}*{h\left\lbrack {n;\theta_{i,h}} \right\rbrack}}} \right)} + {\left( {{c_{i,1}{r^{old}\lbrack n\rbrack}} + {{l^{old}\lbrack n\rbrack}*{h\left\lbrack {n;\theta_{i,h}} \right\rbrack}}} \right) \times \left( {\frac{\left. {c_{i,1}{r\lbrack n\rbrack}} \right)^{2}}{c_{i,1}{r^{old}\lbrack n\rbrack}} + {\sum\limits_{k = 0}^{K - 1}\;{{h\left\lbrack {{n - k};\theta_{i,h}} \right\rbrack}\frac{\left( {l\lbrack k\rbrack} \right)^{2}}{l^{old}\lbrack k\rbrack}}}} \right)}} \right\}}}} & (38)\end{matrix}$It should be observed that, as expected, q_(L)(r^(old), l^(old), c, k₁,k₂, k₃; r^(old), l^(old))=L(r^(old), l^(old), c, k₁, k₂, k₃).

As discussed above, MM algorithms are attractive because they have theproperty that their respective objective functions are monotonicallydecreased with increasing iterations. Thus, it follows from thisproperty that Step 1a can be obtained by solving the followingoptimization problem:

$\begin{matrix}{\left( {r^{({m + 1})},l^{({m + 1})}} \right)\overset{\Delta}{=}{\arg\mspace{14mu}{\min\limits_{{r \geq 0},{l \geq 0}}{q_{L}\left( {r,l,c^{(m)},k_{1}^{(m)},k_{2}^{(m)},{k_{3}^{(m)};r^{(m)}},l^{(m)}} \right)}}}} & (39)\end{matrix}$Before solving equation (39), we make the following definition{circumflex over (d)} _(i) [n;θ _(i,d) ]

c _(i,1) r[n]+l[n]*h[n;θ _(i,h)]  (40)where θ_(i,d)

[r, l, c_(i,1), c_(i,2), k_(i,1), k_(i,2), k_(i,3)]. Now, to determinean update for r and l, we take the partial derivative of the objectivefunction on the right hand side of equation (39) with respect to r andl, respectively, and set the results to zero. For s=0, 1, . . . , N−1the desired partial derivatives are

$\begin{matrix}{\frac{\partial{q_{L}\left( {r,l,c^{(m)},k_{1}^{(m)},k_{2}^{(m)},{k_{3}^{(m)};r^{(m)}},l^{(m)}} \right)}}{\partial{r\lbrack s\rbrack}} = {\sum\limits_{i = 1}^{I}\;\left( {{{- 2}{d_{i}\lbrack n\rbrack}c_{i,1}^{(m)}} + {2{{\hat{d}}_{i}\left\lbrack {n;\theta_{i,d}^{(m)}} \right\rbrack}c_{i,1}^{(m)}\frac{r\lbrack s\rbrack}{r^{(m)}\lbrack s\rbrack}}} \right)}} & (41) \\{\frac{\partial{q_{L}\left( {r,l,c^{(m)},k_{1}^{(m)},k_{2}^{(m)},{k_{3}^{(m)};r^{(m)}},l^{(m)}} \right)}}{\partial{l\lbrack s\rbrack}} = {\sum\limits_{i = 1}^{I}{\sum\limits_{n = 0}^{N - 1}\;\left( {{{- 2}{d_{i}\lbrack n\rbrack}{h\left\lbrack {{n - s};\theta_{i,h}^{(m)}} \right\rbrack}} + {2{{\hat{d}}_{i}\left\lbrack {n;\theta_{i,d}^{(m)}} \right\rbrack}{h\left\lbrack {{n - s};\theta_{i,h}^{(m)}} \right\rbrack}\frac{l\lbrack s\rbrack}{l^{(m)}\lbrack s\rbrack}}} \right)}}} & (42) \\{\mspace{425mu}{= {\sum\limits_{i = 1}^{I}\;\left( {{{- 2}{\sum\limits_{n = 0}^{N - 1}\;{{d_{i}\lbrack n\rbrack}{h\left\lbrack {{n - s};\theta_{i,h}^{(m)}} \right\rbrack}}}} + {2\frac{l\lbrack s\rbrack}{l^{(m)}\lbrack s\rbrack}{\sum\limits_{n = 0}^{N - 1}\;{{{\hat{d}}_{i}\left\lbrack {n;\theta_{i,d}^{(m)}} \right\rbrack}{h\left\lbrack {{n - s};\theta_{i,h}^{(m)}} \right\rbrack}}}}} \right)}}} & (43) \\{\mspace{419mu}{= {\sum\limits_{i = 1}^{I}\;\left\{ {{{- 2}\;{d_{i}\lbrack s\rbrack}*{h\left\lbrack {{- s};\theta_{i,h}^{(m)}} \right\rbrack}} + {2\frac{l\lbrack s\rbrack}{l^{(m)}\lbrack s\rbrack}\left( \;{{{\hat{d}}_{i}\left\lbrack {s;\theta_{i,d}^{(m)}} \right\rbrack}*{h\left\lbrack {{- s};\theta_{i,h}^{(m)}} \right\rbrack}} \right)}} \right\}}}} & (44)\end{matrix}$where θ_(i,d) ^((m))

[r^([m]), l^([m]), c_(i,1) ^((m)), c_(i,2) ^((m)), k_(i,1) ^((m)),k_(i,2) ^((m)), k_(i,3) ^((m))] and θ_(i,h) ^((m))

[c_(i,2) ^((m)), k_(i,1) ^((m)), k_(i,2) ^((m)), k_(i,3) ^((m))], and{circumflex over (d)}_(i)[n; θ_(i,d) ^((m))] is the current estimate ofthe data point d_(i)[n]. Setting equations (41) and (44) to zero leadsto the desired updates for the sampled RV and LV tissue curves:

$\begin{matrix}{{{r^{({m + 1})}\lbrack n\rbrack} = {{r^{(m)}\lbrack n\rbrack}\frac{\sum_{i = 1}^{I}{c_{i,1}^{(m)}{d_{i}\lbrack n\rbrack}}}{\sum_{i = 1}^{I}{c_{i,1}^{(m)}{{\hat{d}}_{i}\left\lbrack {n;\theta_{i,d}^{(m)}} \right\rbrack}}}}},{n = 0},1,\ldots,{N - 1}} & (45) \\{{{l^{({m + 1})}\lbrack n\rbrack} = {{l^{(m)}\lbrack n\rbrack}\frac{\sum_{i = 1}^{I}{{d_{i}\lbrack n\rbrack}*{h\left\lbrack {{- n};\theta_{i,h}^{(m)}} \right\rbrack}}}{\sum_{i = 1}^{I}{{{\hat{d}}_{i}\left\lbrack {n;\theta_{i,d}^{(m)}} \right\rbrack}*{h\left\lbrack {{- n};\theta_{i,h}^{(m)}} \right\rbrack}}}}},{n = 0},1,\ldots,{N - 1}} & (46)\end{matrix}$

The iterates defined by equations (45) and (46) have two importantproperties: (i) they satisfy the desired monotonicity properties givenin equations (35), and (ii) are nonnegative provided the initialestimates, r^((m)), l^((m)), c^((m)), k_(i,1) ^((m)), k_(i,2) ^((m)),and k_(i,3) ^((m)) are nonnegative. In preliminary simulation studies,we have found that the iteration in equation (45) has not consistentlyresulted in accurate estimates for the sampled RV tissue curve. Webelieve the inaccuracy may be due to the fact that the myocardium TACscontains limited information about the RV tissue curve. Therefore,although the update in equation (45) is theoretically meaningful, it maybe advisable to use the estimate for the RV tissue curve obtained fromthe algorithm described below based on a semi-parametric model.

Solution to Step 2:

The minimization in Step 2 is equivalent to the following I minimizationproblems:

     For  i = 1, 2, …  , I,                                           (47)$\begin{matrix}{k_{i,2}^{({m + 1})}\overset{\Delta}{=}{\arg\;{\min\limits_{k_{i,2} \geq 0}{\overset{N - 1}{\sum\limits_{n = 0}}\left\{ {{d_{i}\lbrack n\rbrack} - \left( {{c_{i,1}^{(m)}{r^{({m + 1})}\lbrack n\rbrack}} + {c_{i,2}^{(m)}{l^{({m + 1})}\lbrack n\rbrack}} +} \right.} \right.}}}} \\{{{k_{i,1}^{(m)}\left( {{2/T_{S}} + \left( {k_{i,2} + k_{i,3}^{(m)}} \right)} \right)}^{- 1}{l^{({m + 1})}\lbrack n\rbrack}*{h_{F}\left\lbrack {n;p_{i,2}^{(m)}} \right\rbrack}} +} \\\left. \left. {k_{i,1}^{(m)}{k_{i,3}^{(m)}\left( {2/T_{S}} \right)}^{- 1}\left( {{2/T_{S}} + \left( {k_{i,2} + k_{i,3}^{(m)}} \right)} \right)^{- 1}{l^{({m + 1})}\lbrack n\rbrack}*{h_{T}\left\lbrack {n;p_{i,2}^{(m)}} \right\rbrack}} \right) \right\}^{2}\end{matrix}$                                            (48)$\begin{matrix}{k_{i,3}^{({m + 1})}\overset{\Delta}{=}{\arg\;{\min\limits_{k_{i,2} \geq 0}{\overset{N - 1}{\sum\limits_{n = 0}}\left\{ {{d_{i}\lbrack n\rbrack} - \left( {{c_{i,1}^{({m + 1})}{r^{({m + 1})}\lbrack n\rbrack}} + {c_{i,2}^{(m)}{l^{({m + 1})}\lbrack n\rbrack}} +} \right.} \right.}}}} \\{{{k_{i,1}^{(m)}\left( {{2/T_{S}} + \left( {k_{i,2}^{({m + 1})} + k_{i,3}} \right)} \right)}^{- 1}{l^{({m + 1})}\lbrack n\rbrack}*{h_{F}\left\lbrack {n;p_{i,3}^{({m + 1})}} \right\rbrack}} +} \\\left. \left. {k_{i,1}^{(m)}{k_{i,3}\left( {2/T_{S}} \right)}^{- 1}\left( {{2/T_{S}} + \left( {k_{i,2}^{({m + 1})} + k_{i,3}} \right)} \right)^{- 1}{l^{({m + 1})}\lbrack n\rbrack}*{h_{T}\left\lbrack {n;p_{i,2}^{({m + 1})}} \right\rbrack}} \right) \right\}^{2}\end{matrix}$ where$\mspace{79mu}{p_{i,2}^{(m)}\overset{\Delta}{=}{\left( {{2/T_{S}} + \left( {k_{i,2} + k_{i,3}^{(m)}} \right)} \right)\left( {{2/T_{S}} + \left( {k_{i,2} + k_{i,3}^{(m)}} \right)} \right)^{- 1}\mspace{140mu}(49)}}$$\mspace{79mu}{p_{i,3}^{(m)}\overset{\Delta}{=}{\left( {{2/T_{S}} + \left( {k_{i,2}^{m} + k_{i,3}} \right)} \right)\left( {{2/T_{S}} + \left( {k_{i,2}^{(m)} + k_{i,3}} \right)} \right)^{- 1}\mspace{155mu}(50)}}$Each of the above one-dimensional optimization problems can be solvedusing a line search algorithm such as the golden section method [21].

Solution to Step 3a:

The minimization in Step 3 is equivalent to the following Ione-dimensional minimization problems:

     For  i = 1, 2, …, I,                                          (51)$\begin{matrix}{\theta_{i}^{({m + 1})}\overset{\Delta}{=}{\arg\mspace{14mu}{\min\limits_{\theta_{i} \geq 0}{\overset{N - 1}{\sum\limits_{n = 0}}\left\{ {{d_{i}\lbrack n\rbrack} - \left( {{c_{i,1}{r^{({m + 1})}\lbrack n\rbrack}} + {c_{i,2}{l^{({m + 1})}\lbrack n\rbrack}} +} \right.} \right.}}}} \\{{{k_{i,1}\left( {{2\text{/}T_{s}} + {\overset{\_}{k}}_{i}^{({m + 1})}} \right)}^{- 1}{l^{({m + 1})}\lbrack n\rbrack}*{h_{F}\left\lbrack {n;p_{i}^{({m + 1})}} \right\rbrack}} +} \\\left. \left. {{k_{i,3}^{({m + 1})}\left( {2\text{/}T_{s}} \right)}^{- 1}\left( {{2\text{/}T_{s}} + k_{i}^{({m + 1})}} \right)^{- 1}{l^{({m + 1})}\lbrack n\rbrack}*{h_{T}\left\lbrack {n;p_{i}^{({m + 1})}} \right\rbrack}} \right) \right\}^{2}\end{matrix}$$\mspace{79mu}{{{{where}\mspace{14mu}\theta_{i}}\overset{\Delta}{=}{{\left\lbrack {c_{i,1},c_{i,2},k_{i,1}} \right\rbrack\mspace{14mu}{and}\mspace{14mu}\theta_{i}^{(m)}}\overset{\Delta}{=}{\left\lbrack {c_{i,1}^{(m)},c_{i,2}^{(m)},k_{i,1}^{(m)}} \right\rbrack.\mspace{14mu}{Also}}}},\mspace{79mu}{p_{i}^{(m)}\overset{\Delta}{=}{{\left( {{2\text{/}T_{s}} - k_{i}^{(m)}} \right)\left( {{2\text{/}T_{s}} + k_{i}^{(m)}} \right)^{- 1}\mspace{14mu}{where}\mspace{14mu}{\overset{\_}{k}}_{i}^{(m)}}\overset{\Delta}{=}{k_{i,2}^{(m)} + {{k_{i,3}^{(m)}.\mspace{79mu}{For}}\mspace{20mu}{all}\mspace{14mu} n}}}},{i\mspace{14mu}{let}}}$                                          (52)$\mspace{79mu}{{s_{i,1}^{(m)}\lbrack n\rbrack}\overset{\Delta}{=}{r^{(m)}\lbrack n\rbrack}}$                                          (53)$\mspace{79mu}{{s_{i,2}^{(m)}\lbrack n\rbrack}\overset{\Delta}{=}{l^{(m)}\lbrack n\rbrack}}$                                          (54)${s_{i,3}^{(m)}\lbrack n\rbrack}\overset{\Delta}{=}{{\left( {{2\text{/}T_{s}} + {\overset{\_}{k}}_{i}^{(m)}} \right)^{- 1}{l^{(m)}\lbrack n\rbrack}*{h_{F}\left\lbrack {n;p_{i}^{(m)}} \right\rbrack}} + {{k_{i,3}^{(m)}\left( {2\text{/}T_{s}} \right)}^{- 1}\left( {{2\text{/}T_{s}} + {\overset{\_}{k}}_{i}^{(m)}} \right)^{- 1}{l^{(m)}\lbrack n\rbrack}*{h_{T}\left\lbrack {n;p_{i}^{(m)}} \right\rbrack}}}$It follows that equation (51) can be conveniently expressed as:

$\begin{matrix}\begin{matrix}{\theta_{i}^{({m + 1})}\overset{\Delta}{=}{\arg\mspace{14mu}{\min\limits_{\theta_{i} \geq 0}{\sum\limits_{n = 0}^{N - 1}\;\left\{ {{d_{i}\lbrack n\rbrack} - \left( {{c_{i,1}{s_{i,1}^{({m + 1})}\lbrack n\rbrack}} +} \right.} \right.}}}} \\\left. {{c_{i,2}{s_{i,2}^{({m + 1})}\lbrack n\rbrack}} + {k_{i,1}{s_{i,3}^{({m + 1})}\lbrack n\rbrack}}} \right\}^{2}\end{matrix} & (55) \\{\mspace{65mu}{{= \left. {\arg\mspace{14mu}\min\limits_{\theta_{i} \geq 0}}||{d_{i} - {A_{i}^{({m + 1})}\theta_{i}}}\mathop{\text{||}}_{2}^{2}\mspace{79mu}{where}\mspace{14mu}||{\cdot \left. ||{}_{2}\mspace{14mu}{{is}\mspace{14mu}{the}\mspace{14mu} L_{2}\text{-}{norm}} \right.} \right.},\mspace{79mu}{d_{i} = \left\lbrack {{d_{i}\lbrack 0\rbrack},{d_{i}\lbrack 1\rbrack},\ldots,{d_{i}\left\lbrack {N - 1} \right\rbrack}} \right\rbrack^{T}},{and}}} & (56) \\{\mspace{79mu}{A_{i}^{({m + 1})} = \begin{pmatrix}{s_{i,1}^{({m + 1})}\lbrack 0\rbrack} & {s_{i,2}^{({m + 1})}\lbrack 0\rbrack} & \ldots & {s_{i,3}^{({m + 1})}\lbrack 0\rbrack} \\{s_{i,1}^{({m + 1})}\lbrack 1\rbrack} & {s_{i,2}^{({m + 1})}\lbrack 1\rbrack} & \ldots & {s_{i,3}^{({m + 1})}\lbrack 1\rbrack} \\\vdots & \vdots & \ddots & \vdots \\{s_{i,1}^{({m + 1})}\left\lbrack {N - 1} \right\rbrack} & {s_{i,2}^{({m + 1})}\left\lbrack {N - 1} \right\rbrack} & \ldots & {s_{i,3}^{({m + 1})}\left\lbrack {N - 1} \right\rbrack}\end{pmatrix}}} & (57)\end{matrix}$

Observe that Step 3 requires the solution to equation (56) for alli=1,2, . . . , I. However, since we propose to solve Step 3a instead ofStep 3, we use a majorizing function for a certain class of linear LSobjective functions that was put forth by De Pierro in “On the relationbetween the ISRA and the EM algorithms for positron emissiontomography,” IEEE Transactions on Medical Imaging, vol. 12, pp. 328-333,1993, which is incorporated by reference. Given that θ_(i), d_(i), andthe matrix A_(i) ^((m+1)) are nonnegative, De Pierro's result can beused to obtain the following majorizing function for ∥d_(i)−A_(i)^((m+1))θ_(i)∥₂ ² about the point θ_(i) ^((m))

$\begin{matrix}{{q_{1}\left( {\theta_{i},\theta_{i}^{(m)}} \right)}\overset{\Delta}{=}{{\sum\limits_{j = 0}^{N - 1}\;{d_{i}\lbrack j\rbrack}^{2}} - {2{\sum\limits_{j = 0}^{N - 1}\;{{d_{i}\lbrack j\rbrack}\left\lbrack {A_{i}^{({m + 1})}\theta_{i}} \right\rbrack}_{j}}} + {\sum\limits_{j = 0}^{N - 1}\;{\sum\limits_{k = 1}^{3}\;{\left\lbrack {A_{i}^{({m + 1})}\theta_{i}^{(m)}} \right\rbrack_{j}A_{i,j,k}^{({m + 1})}\frac{\theta_{i,k}^{2}}{\theta_{i,k}^{(m)}}}}}}} & (58)\end{matrix}$where A_(i,j,k) ^((m+1)) is the element in j^(th) row and k^(th) columnof the matrix A_(i) ^((m+1)), and θ_(i,k) is the k^(th) element of thevector θ_(i). Additionally, for a vector ν the quantity [ν]k is definedto be the k^(th) element of ν. Now, the partial derivative of themajorizing function q₁ with respect to θ_(i,k) equals

$\begin{matrix}{\frac{\partial{q_{1}\left( {\theta_{i},\theta_{i}^{(m)}} \right)}}{\partial\theta_{i,k}} = {{- {2\left\lbrack {\left( A_{i}^{({m + 1})} \right)^{T}d_{i}} \right\rbrack}_{k}} + {2{\frac{\theta_{i,k}}{\theta_{i,k}^{(m)}}\left\lbrack {\left( A_{i}^{({m + 1})} \right)^{T}\left( {A_{i}^{({m + 1})}\theta_{i}^{(m)}} \right)} \right\rbrack}_{k}}}} & (59)\end{matrix}$

Setting this partial derivative to zero yields the desired update, whichis also known as the iterative space image reconstruction algorithm(IRSA) as described in M. E. Daube-Witherspoon and G. Muehllehner,“Treatment of axial data in three-dimensional PET,” Journal of NuclearMedicine, vol. 28, pp. 1717-1724, 1987, which is incorporated byreference.

$\begin{matrix}{{\theta_{i,k}^{({m + 1})} = {{\theta_{i,k}^{(m)}\frac{\left\lbrack {\left( A_{i}^{({m + 1})} \right)^{T}d_{i}} \right\rbrack_{k}}{\left\lbrack {\left( A_{i}^{({m + 1})} \right)^{T}\left( {A_{i}^{({m + 1})}\theta_{i}^{(m)}} \right)} \right\rbrack_{k}}\mspace{14mu} k} = 1}},2,3} & (60)\end{matrix}$

From equation (60), the updates for the parameters c_(i,1), c_(i,2), andk_(i,1) can be equivalently written as

$\begin{matrix}{{c_{i,k}^{({m + 1})} = {{c_{i,k}^{(m)}\frac{\left\lbrack {\left( A_{i}^{({m + 1})} \right)^{T}d_{i}} \right\rbrack_{k}}{\left\lbrack {\left( A_{i}^{({m + 1})} \right)^{T}\left( {A_{i}^{({m + 1})}\theta_{i}^{(m)}} \right)} \right\rbrack_{k}}\mspace{14mu} k} = 1}},2} & (61) \\{k_{i,1}^{({m + 1})} = {k_{i,1}^{(m)}\frac{\left\lbrack {\left( A_{i}^{({m + 1})} \right)^{T}d_{i}} \right\rbrack_{3}}{\left\lbrack {\left( A_{i}^{({m + 1})} \right)^{T}\left( {A_{i}^{({m + 1})}\theta_{i}^{(m)}} \right)} \right\rbrack_{3}}}} & (62)\end{matrix}$

C. Penalized MBF Algorithm

When the myocardium TACs are noisy, the MBF estimates generated by theiteration in equation (62) may be so noisy that some smoothing isnecessary. Our approach for the addressing the noise is to add alimiting factor or penalty function that penalizes solutions where theMBF estimates for adjacent voxels differ significantly. Specifically, wedevelop the following penalized least-squares (PLS) method:

$\begin{matrix}{\theta_{i}^{({m + 1})}\overset{\Delta}{=}\left. {\arg\mspace{14mu}\min\limits_{\theta_{i} \geq 0}}||{d_{i} - {A_{i}^{({m + 1})}\theta_{i}}}\mathop{\text{||}}_{2}^{2}{{+ \lambda}\;{P\left( k_{1} \right)}} \right.} & (63)\end{matrix}$where P(k₁) is a penalty function, N_(i) is a neighborhood about thei^(th) voxel, and the penalty parameter, λ, determines the level ofinfluence of the penalty function. For the penalty function, we choosethe following quadratic penalty function:

$\begin{matrix}{{P\left( k_{1} \right)} = {\sum\limits_{i = 1}^{I}\;{\sum\limits_{j \in N_{i}}\left( {k_{i,1} - k_{j,1}} \right)^{2}}}} & (64)\end{matrix}$

A choice for N_(i) is the intersection of the set of eight nearestvoxels to the i^(th) voxel and the set of voxels that lie in themyocardium. Note, penalty functions could also be used to enforcesmoothness on the other parameters c, k₂, and k₃.

The following majorizing function for the penalty function P about thepoint k_(i,1)(m) was provided in G. Wang and J. Qi, “Penalizedlikelihood PET image reconstruction using patch-based edge-preservingregularization,” IEEE Transactions on Medical Imaging, vol. 31, pp.2194-2204, 2012, which is incorporated herein by reference:

$\begin{matrix}{{q_{2}\left( {k_{1},k_{1}^{(m)}} \right)}\overset{\Delta}{=}{\sum\limits_{i = 1}^{I}\;{\sum\limits_{j \in N_{i}}\left\lbrack {{2\left( {k_{i,1} - \frac{k_{i,1}^{(m)} + k_{j,1}^{(m)}}{2}} \right)^{2}} + {2\left( {k_{i,j} - \frac{k_{i,1}^{(m)} + k_{j,1}^{(m)}}{2}} \right)^{2}}} \right\rbrack}}} & (65)\end{matrix}$

Because q₁ is a majorizing function for the LS objective function∥d_(i)−A_(i) ^((m+1)) θ_(i)∥₂ ², it follows thatq(θ_(i),θ_(i) ^((m)))

q ₁(θ_(i),θ_(i) ^((m)))+λq ₂(k ₁ ,k ₁ ^((m)))  (66)is a majorizing function for the PLS objective function in equation (63)about the point θ_(i) ^((m)). Consequently, an MM algorithm forobtaining PLS estimates of the MBF parameters can be expressed as

$\begin{matrix}{\theta_{i}^{({m + 1})}\overset{\bigtriangleup}{=}{\arg\;{\min\limits_{\theta_{i} \geq 0}{q\left( {\theta_{i},\theta_{i}^{(m)}} \right)}}}} & (67)\end{matrix}$

To get the update for k_(i,1) using the PLS formulation, we first takethe derivative of q with respect to k_(i,1):

$\begin{matrix}{\frac{\partial{q\left( {\theta_{i},\theta_{i}^{(m)}} \right)}}{\partial k_{i,1}} = {{- {2\left\lbrack {\left( A_{i}^{({m + 1})} \right)^{T}d_{i}} \right\rbrack}_{3}} + {2{\frac{k_{i,1}}{k_{i,1}^{(m)}}\left\lbrack {\left( A_{i}^{({m + 1})} \right)^{T}\left( {A_{i}^{({m + 1})}\theta_{i}^{(m)}} \right)} \right\rbrack}_{3}} + {8\lambda\; k_{i,1}{\quad\left| N_{i} \middle| {{- 8}\lambda{\sum\limits_{j \in N_{i}}\left\lbrack \left( \frac{k_{i,1}^{(m)} + k_{j,1}^{(m)}}{2} \right) \right\rbrack}} \right.}}}} & (68)\end{matrix}$where |N_(i)| is the number of elements in the set N_(i). Now, equatingthis result to zero leads to the desired result:

$\begin{matrix}{{k_{i,1}^{({m + 1})} = \frac{{k_{i,1}^{(m)}\left\lbrack {\left( A_{i}^{({m + 1})} \right)^{T}d_{i}} \right\rbrack}_{3} + {2\lambda\; k_{i,1}^{(m)}{\Sigma_{j \in N_{i}}\left( {k_{i,1}^{(m)} + k_{j,1}^{(m)}} \right)}}}{\left. {\left\lbrack {\left( A_{i}^{({m + 1})} \right)^{T}\left( {A_{i}^{({m + 1})}\theta_{i}^{(m)}} \right)} \right\rbrack_{3} + {4\lambda\; k_{i,1}^{(m)}}} \middle| N_{i} \right|}},{i = 1},2,\ldots,I} & (69)\end{matrix}$The penalty function does not depend on the parameters c_(i,1) andc_(i,2) so the updates c_(i,1) ^((m+1)) and c_(i,1) ^((m+1)) are givenby equation (61).

In the PLS method, a choice is made for the penalty parameter λ. In ourexperiments, we chose the penalty parameter experimentally. However, thepopular L-curve method could be used instead to determine the penaltyparameter as discussed in P. C. Hansen, “Analysis of discrete ill-posedproblems by means of the L-curve,” SIAM Review, vol. 34, pp. 561-580,1992; P. C. Hansen and D. P. OLeary, “The use of the L-curve in theregularization of discrete ill-posed problems,” SIAM Journal ScientificComputing, vol. 14, p. 14871503, 1993; and T. Reginska, “Aregularization parameter in discrete ill-posed problems,” SIAM JournalScientific Computing, vol. 17, p. 740749, 1996, each of which areincorporated herein in their entireties.

1. Semi-Parametric Model for Sampled RV and LV Tissue Curves

From an understanding of human physiology, it is well known that the RVand LV tissue curves will decay to zero. With the expectation of moreaccurate MBF estimation, we incorporate this knowledge by modeling theLV tissue curve as

$\begin{matrix}{{s_{l}\left\lbrack \;{{n;n_{0,l}},\beta} \right\rbrack}\overset{\Delta}{=}\left\{ \begin{matrix}{{l\lbrack n\rbrack},} & {0 \leq n \leq n_{0,l}} \\{{{l\left\lbrack n_{o,l} \right\rbrack}\beta^{({n - n_{0,l}})}},} & {n_{o,i} < n \leq {N - 1}} \\{0,} & {{otherwise},}\end{matrix} \right.} & (70)\end{matrix}$where n_(0,l) and β are unknown constants. We refer to this model as asemi-parametric model because only the sampled LV tissue curve valuesfor n≥n_(0,l) are described by a parametric model. In oursemi-parametric model based method, we account for the fact that themyocardium TACs contain little information about the RV tissue curve bycombining the initial RV tissue curve estimate, r(0), which is estimatedin a non-parametric fashion, with the semi-parametric model for the LVtissue curve in the following manner:

$\begin{matrix}{{s_{r}\left\lbrack {{n;n_{0,r}},\beta} \right\rbrack}\overset{\Delta}{=}\left\{ \begin{matrix}{{r^{(0)}\lbrack n\rbrack},} & {0 \leq n \leq n_{0,r}} \\{{{r^{(0)}\left\lbrack n_{0,r} \right\rbrack}\beta^{({n - n_{0,r}})}},} & {n_{0,r} < n \leq {N - 1}} \\{0,} & {{otherwise},}\end{matrix} \right.} & (71)\end{matrix}$where the parameter n_(0,r) is unknown. For the estimated sampled RV andLV tissue curves to be nonnegative and decay to zero, the parameter βmust satisfy the constraint 0<β<1.

2. LS Estimation of Sampled RV and LV Tissue Curves for n_(0,r) andn_(0,l) Known

For the moment, we will assume that the parameters n_(0,r) and n_(0,l)are known. Using the models for the sampled RV and LV tissue curvesexpressed by equations (70) and (71), the associated LS objectivefunction is given by

$\begin{matrix}{{\overset{\_}{L}\left( {l_{n_{0,l}},\beta,c,k_{1},k_{2},k_{3}} \right)} = {\sum\limits_{i = 1}^{I}\;{\sum\limits_{n = 0}^{N - 1}\;\left\{ {{d_{i}\lbrack n\rbrack} - \left( {{c_{i,1}{s_{r}\left\lbrack {{n;n_{0,r}},\beta} \right\rbrack}} + {{s_{l}\left\lbrack {{n;n_{0,l}},\beta} \right\rbrack}*{h\left\lbrack {n;\theta_{i,h}} \right\rbrack}}} \right)} \right\}^{2}}}} & (72)\end{matrix}$where l_(n) _(o,l)

[l[0], l[1], . . . , l[n_(o,l)]]. Therefore, an alternative approach toStep 1a are the following two steps:

Model Based Step 1a: Find an iterate l_(n) _(0,l) ^((m+1)) such that:L (l _(n) _(0,l) ^((m+1)),β^((m)) ,c ^((m)) ,k ₁ ^((m)) ,k ₂ ^((m)) ,k ₃^((m)))≤ L (l _(n) _(0,l) ^((m)),β^((m)) ,c ^((m)) ,k ₁ ^((m)) ,k ₂^((m)) ,k ₃ ^((m)))  (73)

Model Based Step 1b

$\begin{matrix}{\beta^{({m + 1})}\overset{\Delta}{=}{\arg\mspace{14mu}{\min\limits_{\epsilon \leq \beta \leq {1 - \epsilon}}{\overset{\_}{L}\left( {l_{n_{0,l}}^{({m + 1})},\beta^{(m)},c^{(m)},k_{1}^{(m)},k_{2}^{(m)},k_{3}^{(m)}} \right)}}}} & (74)\end{matrix}$where ∈ is a small positive constant. Note replacing the constraint0<β<1 with ∈≤β≤1−∈ makes the problem in Model Based Step 1b welldefined.

First, we address the problem in Model Based Step 1a. Let the N×1vectors s_(r) ^(old) and s_(l) ^(old) denote an estimate of thesemi-parametric RV and LV tissue curves, respectively:s _(r) ^(old)

[s _(r) ^(old)[0;n _(0,r),β^(old) ],s _(r) ^(old)[1;n _(0,r),β^(old) ],. . . ,s _(r) ^(old) [N−1;n _(0,r),β^(old)]]  (75)=[r ⁽⁰⁾[0],r ⁽⁰⁾[1], . . . ,r ⁽⁰⁾ [n _(0,r) ],r ⁽⁰⁾ [n _(0,r)](β^(old))¹, . . . ,r ⁽⁰⁾ [n _(0,r)](β^(old))^(N−1−n) ^(0,r) ]  (76)s _(l) ^(old)

[s _(l) ^(old)[0;n _(0,l),β^(old) ],s _(l) ^(old)[1;n _(0,l),β^(old) ],. . . ,s _(l) ^(old) [N−1;n _(0,l),β^(old)]]  (77)=[l ^(old)[0],l ^(old)[1], . . . ,l ^(old) [n _(0,l) ],l ^(old) [n_(0,l)](β^(old))¹ , . . . ,l ^(old) [n _(0,l)](β^(old))^(N−1−n) ^(0,l)]  (78)Mimicking the steps used to derive equation (38), a majorizing functionfor L at (s_(r) ^(old), s_(l) ^(old)) is:

$\begin{matrix}{{q_{\overset{\_}{L}}\left( {{l_{n_{0,l}}\beta},c,k_{1},k_{2},{k_{3};s_{r}^{old}},s_{l}^{old}} \right)} = {\sum\limits_{i = 1}^{I}\;{\sum\limits_{n = 0}^{N - 1}\;\left\{ {{d_{i}\lbrack n\rbrack}^{2} - {2{d_{i}\lbrack n\rbrack}\left( {{c_{i,1}{s_{r}\left\lbrack {{n;n_{0,r}},\beta} \right\rbrack}} + {{s_{l}\left\lbrack {{n;n_{0,l}},\beta} \right\rbrack}*{h\left\lbrack {n;\theta_{i,h}} \right\rbrack}}} \right)} + {\left( {{c_{i,1}s_{r}^{old}\left\{ {{n;n_{0,r}},\beta^{old}} \right\rbrack} + {{s_{l}^{old}\left\lbrack {{n;n_{0,l}},\beta^{old}} \right\rbrack}*{h\left\lbrack {n;\theta_{i,h}} \right\rbrack}}} \right) \times \left( {\frac{\left( {c_{i,1}{s_{r}\left\lbrack {{n;n_{0,r}},\beta} \right\rbrack}} \right)^{2}}{c_{i,1}{s_{r}^{old}\left\lbrack {{n;n_{0,r}},\beta^{old}} \right\rbrack}} + {\sum\limits_{k = 0}^{N - 1}\;{{h\left\lbrack {{n - k};\theta_{i,h}} \right\rbrack}\frac{\left( {s_{l}\left\lbrack {{k;n_{0,l}},\beta} \right\rbrack} \right)^{2}}{s_{l}^{old}\left\lbrack {{k;n_{0,l}},\beta^{old}} \right\rbrack}}}} \right)}} \right\}}}} & (79)\end{matrix}$

For ν=0, 1, . . . , n_(0,l), the derivative of q _(L) with respect tos_(l)[ν; n_(0,l), β_(l)] follows from equation (44) and is equal to

$\begin{matrix}{\frac{\partial{q_{\overset{\_}{L}}\left( {l_{n_{0,l}}^{(m)},\beta^{(m)},c^{(m)},k_{1}^{(m)},k_{2}^{(m)},{k_{3}^{(m)};s_{r}^{(m)}},s_{l}^{(m)}} \right)}}{\partial{s_{l}\left\lbrack {{v;n_{0,l}},\beta} \right\rbrack}} = {\sum\limits_{i = 1}^{I}\;{\sum\limits_{n = 0}^{N - 1}\;\left( {{{- 2}{d_{i}\lbrack n\rbrack}{h\left\lbrack {{n - v};\theta_{i,h}^{(m)}} \right\rbrack}} + {2{{\overset{\_}{d}}_{i}\left\lbrack {{n;s_{r}^{(m)}},s_{l}^{(m)},c_{i,1}^{(m)},\theta_{i,h}^{(m)}} \right\rbrack}{h\left\lbrack {{n - v};\theta_{i,h}^{(m)}} \right\rbrack}\frac{s_{l}\left\lbrack {{v;n_{0,l}},\beta^{(m)}} \right\rbrack}{s_{l}^{(m)}\left\lbrack {{v;n_{0,l}},\beta^{(m)}} \right\rbrack}}} \right)}}} & (80) \\{= {\sum\limits_{i = 1}^{I}\;\left\{ {{{- 2}{d_{i}\lbrack v\rbrack}*{h\left\lbrack {{- v};\theta_{i,h}^{(m)}} \right\rbrack}} + {2\frac{s_{l}\left\lbrack {v;{n_{0,l}\beta^{(m)}}} \right\rbrack}{s_{l}^{(m)}\left\lbrack {{v;n_{0,l}},\beta^{(m)}} \right\rbrack}\left( {{{\overset{\_}{d}}_{i}\left\lbrack {{v;s_{r}^{(m)}},s_{l}^{(m)},c_{i,1}^{(m)},\theta_{i,h}^{(m)}} \right\rbrack}*{h\left\lbrack {{- v};\theta_{i,h}^{(m)}} \right\rbrack}} \right)}} \right\}}} & (81)\end{matrix}$where s_(r) ^((m)) and s_(l) ^((m)) are derived in a similar way ass_(r) ^((old)) and s_(l) ^((old)), respectively, andd _(i) [n;s _(r) ^((m)) ,s _(l) ^((m)) ,c _(i,1) ^((m)),θ_(i,h) ^((m)) ]

c _(i,1) ^((m)) s _(r) ^((m)) [n;n _(0,r),β^((m)) ]+s _(l) ^((m)) [n;n_(0,l),β^((m)) ]*h[n;θ _(i,h) ^((m))]  (82)is an estimate for the data point d_(i)[n] using the semi-parametricmodel for the sampled RV and LV tissue curves in equations (70) and(71). From the definition in equation (70), for n=0, 1 . . . . ,n_(0,l), it follows that s_(l)[ν; n_(0,l), β]=l[ν] and s_(l) ^((m))[ν;n_(0,l), β^((m))]=l^((m))[ν]. Thus, for n=0, 1, . . . , n_(0,l), thepartial derivative in equation (82) becomes

$\begin{matrix}{\frac{\partial{q_{\overset{\_}{L}}\left( {l_{n_{0,l}}^{(m)},\beta^{(m)},c^{(m)},k_{1}^{(m)},k_{2}^{(m)},{k_{3}^{(m)};s_{r}^{(m)}},s_{l}^{(m)}} \right)}}{\partial{s_{l}\left\lbrack {{v;n_{0,l}},\beta} \right\rbrack}} = \frac{\partial{q_{\overset{\_}{L}}\left( {l_{n_{0,l}},\beta^{(m)},c^{(m)},k_{1}^{(m)},k_{2}^{(m)},{k_{3}^{(m)};s_{r}^{(m)}},s_{l}^{(m)}} \right)}}{\partial{l\lbrack v\rbrack}}} & (83) \\{= {\sum\limits_{i = 1}^{I}\;\left\{ {{{- 2}{d_{i}\lbrack v\rbrack}*{h\left\lbrack {{- v};\theta_{i,h}^{(m)}} \right\rbrack}} + {2\frac{l\lbrack v\rbrack}{l^{(m)}\lbrack v\rbrack}\left( {{{\overset{\_}{d}}_{i}\left\lbrack {{v;s_{r}^{(m)}},s_{l}^{(m)},c_{i,l}^{(m)},\theta_{i,h}^{(m)}} \right\rbrack}*{h\left\lbrack {{- v};\theta_{i,h}^{(m)}} \right\rbrack}} \right)}} \right\}}} & (84)\end{matrix}$Setting equation (85) to zero leads to the desired update for thesampled LV tissue curve for n=0, 1, . . . , n_(0,l)

$\begin{matrix}{{l^{({m + 1})}\lbrack n\rbrack} = {{l^{(m)}\lbrack n\rbrack}\frac{\sum_{i = 1}^{I}{{d_{i}\lbrack n\rbrack}*{h\left\lbrack {{- n};\theta_{i,h}^{(m)}} \right\rbrack}}}{\sum_{i = 1}^{I}{{{\overset{\_}{d}}_{i}\left\lbrack {{n;s_{r}^{(m)}},s_{l}^{(m)},c_{i,1}^{(m)},\theta_{i,h}^{(m)}} \right\rbrack}*{h\left\lbrack {{- n};\theta_{i,h}^{(m)}} \right\rbrack}}}}} & (85)\end{matrix}$

The next update for β is obtained by solving the minimization problem inModel Based Step 1b by applying a one-dimensional line search such asdescribed by M. S. Bazaraa, H. D. Sherali, and C. M. Shetty, NonlinearProgramming Theory and Algorithms. New York, N.Y.: John Wiley & Sons,Second Edition, Chapter 8, 1993, which is incorporated by reference.Finally, given β^((m+1)), the next estimates for the RV and LV tissuecurve at the same point n=n_(0,r)+1, n_(0,r)+2, . . . N−1 andn=n_(0,l)+1, n_(0,l)+2, . . . N−1, respectively, are given by:r ^((m+1)) [n]=r ⁽⁰⁾ [n _(0,r)](β^((m+1)))^((n−n) ^(0,r) ⁾  (86)l ^((m+1)) [n]=l ^((m+1)) [n _(0,l)](β^((m+1)))^((n−n) ^(0,l) ⁾  (87)Summarizing, the updates for the sampled RV and LV tissue curves whenthe semi-parametric model is assumed are given by

$\begin{matrix}{{r^{({m + 1})}\lbrack n\rbrack} = \left\{ \begin{matrix}{{r^{(0)}\lbrack n\rbrack},} & {0 \leq n \leq n_{0,r}} \\{{{r^{(0)}\left\lbrack n_{0,r} \right\rbrack}\left( \beta^{({m + 1})} \right)^{({n - n_{0,r}})}},} & {n_{0,r} < n \leq {N - 1}} \\{0,} & {{otherwise},}\end{matrix} \right.} & (88) \\{{l^{(m)}\lbrack n\rbrack} = \left\{ \begin{matrix}{{{l^{(m)}\lbrack n\rbrack}\frac{\sum_{i = 1}^{I}{{d_{i}\lbrack n\rbrack}*{h\left\lbrack {{- n};\theta_{i,h}^{(m)}} \right\rbrack}}}{\sum_{i = 1}^{I}{{{\overset{\_}{d}}_{l}\left\lbrack {{n;s_{r}^{(m)}},s_{l}^{(m)},c_{i,1}^{(m)},\theta_{i,h}^{(m)}} \right\rbrack}*{h\left\lbrack {{- n};\theta_{i,h}^{(m)}} \right\rbrack}}}},} & {0 \leq n \leq n_{0,l}} \\{{{l^{({m + 1})}\left\lbrack n_{o,l} \right\rbrack}\left( \beta^{({m + 1})} \right)^{({n - n_{0,l}})}},} & {n_{0,l} < n \leq {N - 1}} \\{0,} & {otherwise}\end{matrix} \right.} & (89)\end{matrix}$

3. Estimation of the Parameters: n_(0,r) and n_(0,l)

We define d_(avg,r)[n] and d_(avg,l)[n] to be the average of the RV andLV TACs, respectively. Let n_(r,max) and n_(l,max) equal the time pointswhere d_(avg,r)[n] and d_(avg,l)[n] equal their maximum values,respectively. Additionally, let n_(r,half) represent the time pointafter n_(r,max) where d_(avg,r)[n] is halfway down from its maximumvalue (note: n_(l,half) is defined similarly). A reasonable assumptionis that the time points where the sampled RV and LV tissues curves beginto decay exponentially are n_(r,half) and n_(l,half), respectively.

Another approach would be to first compute the K-FADS-II parameters foreach (n_(0,r), n_(0,l)) pair that comes from a set of (n_(ϕ,r), n_(ϕ,l))pairs, such as the set{(n _(r,half) ,n _(l,half)),(n _(r,half)+1,n _(l,half)+1), . . . ,(n_(r,half) +T _(max) ,n _(l,half) +T _(max))}  (90)where T_(max) would be specific by the user. Then, an error criterionknown as the minimum description criterion such as that described in B.Porat, Digital Processing of Random Signals: Theory and Methods. NewYork, N.Y.: Prentice Hall, 1994, which is incorporated by reference,would be used to determine the optimal estimates for n_(0,r), andn_(0,l), and then the K-FADS-II parameters would be estimated using theproposed algorithm.

E. Initialization Procedure

In this section, we first discuss a straightforward way to obtaininitial estimates of the RV and LV tissue curves. Then, given theseinitial estimates, we present a reliable way to generate initialestimates for c, k₁, k₂, and k₃.

1. Initial RV and LV Tissue Curves

Due to the motion of the heart and the finite resolution of PET imaging,the RV and LV TACs are corrupted by activity in the myocardium so theydo not decay to zero. Keeping in mind the reasoning behind thesemi-parametric model for the RV and LV tissue curves, the initial RVand LV tissue curves are chosen to be

$\begin{matrix}{{r^{(0)}\lbrack n\rbrack} = \left\{ \begin{matrix}{{d_{{avg},r}\lbrack n\rbrack},} & {0 \leq n \leq n_{0,r}} \\{{{d_{{avg},r}\left\lbrack n_{0,r} \right\rbrack}\left( \beta^{(0)} \right)^{({n - n_{0,r}})}},} & {n_{0,r} < n \leq {N - 1}} \\{0,} & {{otherwise},}\end{matrix} \right.} & (91) \\{{l^{(0)}\lbrack n\rbrack} = \left\{ \begin{matrix}{{d_{{avg},l}\lbrack n\rbrack},} & {0 \leq n \leq n_{0,l}} \\{{{d_{{avg},l}\left( n_{0,l} \right\rbrack}\left( \beta^{(0)} \right)^{({n - n_{0,l}})}},} & {n_{o,l} < n \leq {N - 1}} \\{0,} & {otherwise}\end{matrix} \right.} & (92)\end{matrix}$where β⁽⁰⁾ is the initial estimate of the parameter β.

An extension of this initialization procedure is to first determine themaximum value of each of RV TAC and then scale r⁽⁰⁾[n] in equation (91)so that its maximum value equals the average of the Q largest maximumvalues. The corresponding scale factor would be applied to l⁽⁰⁾[n] inequation (92). The underlying assumption behind this extension is thatthe RV and LV TACs with the largest maximum values have the least amountof noise due to activity in the myocardium. Note, in the experimentsdiscussed in Section VIII we used β⁽⁰⁾=0.95 and Q=3.

2. Initial Estimates for the Parameters: k₂, and k₃

Now, we will describe a method for finding initial estimates for k₂ andk₃ given the initial estimates r⁽⁰⁾ and l⁽⁰⁾. The model in equation (14)can be written as:

$\begin{matrix}{{D_{i}(z)} = {{c_{i,1}{R(z)}} + {{L(z)}\frac{\begin{matrix}{{{c_{i,2}\left( {1 - {p_{i}z^{- 1}}} \right)}\left( {1 - z^{- 1}} \right)} +} \\{{{{\overset{\_}{c}}_{i,3}\left( {1 - z^{- 1}} \right)}\left( {1 + z^{-}} \right)} + {{\overset{\_}{c}}_{i,4}\left( {1 + z^{- 1}} \right)}^{2}}\end{matrix}}{\left( {1 - {p_{i}z^{- 1}}} \right)\left( {1 - z^{- 1}} \right)}}}} & (93)\end{matrix}$where c _(i,3)

c_(i,3)k_(i,1)(2/T_(s)+k_(i))⁻¹ and c _(i,4)

c_(i,4)k_(i,1)k_(i,3)(2/T_(s)+k_(i))⁻¹(2/T_(s))⁻¹. Thus, for somesecond-order polynomial B_(i)(z)=b_(i)[0]+b_(i)[1]z⁻¹+b_(i)[2]z⁻² andfirst-order polynomial M_(i)(z)=m_(i)[0]+m_(i)[1]z⁻¹, and first orderpolynomial A_(i)(z)=a_(i)[0]+a_(i)[1]z⁻¹ with a_(i)[0]=1, it followsthat:

$\begin{matrix}{{{D_{i}(z)} = \frac{{{M_{i}(z)}{R(z)}} + {{\overset{\_}{L}(z)}{B_{i}(z)}}}{A_{i}(z)}}{where}} & (94) \\{{\overset{\_}{L}(z)} = \frac{L(z)}{1 - z^{- 1}}} & (95)\end{matrix}$Comparing equations (93) and (94) it can be seen that a_(i)[1]=p_(i).Given the myocardium TACs {d_(i)[n]} and initial estimates for thesampled RV and LV tissue curves, the problem is to estimate theparameter vector ϕ

[a_(i)[1], b_(i)[0], b_(i)[1], b_(i)[2], m_(i)[0], m_(i)[1]].

To estimate ϕ, we propose an equation-error method that is amodification of the well-known Steiglitz-McBride algorithm and based onD _(i)(z)A _(i)(z)=M _(i)(z)R(z)+ L (z)B _(i)(z)  (96)which follows from equation (94). The advantage of using equation (95)instead of equation (94) is that the former equation leads to a linearLS method while the latter results in a nonlinear LS formulation. Thecorresponding time-domain expression for equation (96) isd _(i) [n]*a _(i) [n]=r[n]*m _(i) [n]+l[n]*u[n]*b _(i) [n]  (97)where u[n] is the discrete-time unit step function. It follows fromequation (97) that the parameter ϕ could be estimated by minimizingL_(i,init)(ϕ) where

$\begin{matrix}{{L_{i,{init}}(\phi)}\overset{\Delta}{=}{\sum\limits_{n = 0}^{N - 1}\;\left\{ {{{d_{i}\lbrack n\rbrack}*{a_{i}\lbrack n\rbrack}} - \left( {{{\overset{\_}{r}\lbrack n\rbrack}*{m_{i}\lbrack n\rbrack}} + {{\overset{\_}{l}\lbrack n\rbrack}*{b_{i}\lbrack n\rbrack}}} \right)} \right\}^{2}}} & (98)\end{matrix}$where r[n]

r⁽⁰⁾[n] and l[n]

l⁽⁰⁾[n]*u[n]. Instead, we modify a method suggested by Steiglitz andMcBride for estimating the parameters of the rational system function ofa discrete-time, linear time-invariant system given input-output datafrom the system.

Algorithm 1 below is a summary of the method proposed for determininginitial values for the parameters p

[p₁, p₂, . . . , p_(l)], k₂, and k₃. The reader should keep in mind thefollowing relationships from above: k

k_(i,2)+k_(i,3) and p_(i)

(2/T_(S)−k_(i))(2/T_(S)−k_(i))⁻¹.

Algorithm 1 Initialization Procedure for the Parameteres p, k₂, and k₃Get d_(avg,r)[n] and d_(avg,l)[n] for n = 0, 1, ..., N − 1 n_(0,r) ←n_(r,half) n_(0,l) ← n_(l,half) Choose 0 < β⁽⁰⁾ < 1 Get r⁽⁰⁾ and l⁽⁰⁾using equations (92) and (93) For all i = 1, 2, . . ., I do (99)$\left. \phi_{i}^{(0)}\leftarrow{\arg\mspace{11mu}{\min\limits_{\phi_{i}}\;{L_{i,{init}}\left( \phi_{i} \right)}}} \right.,$ For all j = 0, 1, . . ., J₁ - 1 Get the sequence h[n; a^((s))[1]]${{where}\mspace{14mu}{h\left\lbrack {n;\alpha} \right\rbrack}}\overset{\Delta}{=}{{\mathcal{Z}^{- 1}\left\lbrack \frac{1}{1 - {\alpha\; z^{- 1}}} \right\rbrack}\mspace{14mu}\text{=}\text{∝}^{n}\mspace{14mu}{u\lbrack n\rbrack}}$Compute d_(i) ^((s))[n] ← d_(i)[n] * h[n; a_(i) ^((s))[1]] Compute r^((s))[n] ← r[n] * h[n; a_(i) ^((s))[1]] Compute l ^((s))[n] ← l[n] *h[n; a_(i) ^((s))[1]] Get new estimate for φ₁ (100)${\phi_{i}^{({s + 1})}\overset{\Delta}{=}{\arg\mspace{11mu}{\min\limits_{\phi_{i}}\;{L_{i,{nit}}^{(s)}\left( \phi_{i} \right)}}}},$${{where}\mspace{14mu}{L_{i,{init}}^{(s)}\left( \phi_{i} \right)}}\overset{\Delta}{=}{\sum\limits_{n = 0}^{N - 1}\left( {{{d_{i}^{(s)}\lbrack n\rbrack}*{a_{i}\lbrack n\rbrack}} - \left( {{{{\overset{\_}{r}}^{(s)}\lbrack n\rbrack}*{m_{i}\lbrack n\rbrack}} + {{{\overset{\_}{l}}^{(s)}\lbrack n\rbrack}*{b_{i}\lbrack n\rbrack}}} \right)} \right)^{2}}$ end for p_(i) ⁽⁰⁾ ← â[1] where {circumflex over (φ)}_(i) 

 (â_(i)[1], {circumflex over (b)}_(i)[0], {circumflex over (b)}_(i)[1], 

[2], {circumflex over (m)}_(i)[0], {circumflex over (m)}_(i)[1]) is theresulting estimate for φ_(i)$\left. k_{i}^{(0)}\leftarrow{\frac{2}{T_{S}}\frac{1 - p_{i}^{(0)}}{1 + p_{i}^{(0)}}} \right.$$\left. k_{i,2}^{(0)}\leftarrow{k_{i}^{(0)}\frac{factor}{1 + {factor}}} \right.,$$\left. k_{i,3}^{(0)}\leftarrow{k_{i}^{(0)}\frac{1}{1 + {factor}}} \right.,$end for

The optimization problem in (99) and (100) are straightforward becausethey are in the form of an unconstrained, linear LS estimation problem,which has a known solution. Specifically, the problems are of the form:

$\begin{matrix}{\hat{x}\overset{\Delta}{=}\left. {\arg\mspace{14mu}\min\limits_{x}}||{y - {Ax}}||_{2}^{2} \right.} & (101)\end{matrix}$where y is the data and A is a known matrix. The unique minimum L₂ normsolution is given by {circumflex over (x)}=A+b, where A⁺ is thepseudoinverse of A per known applied linear algebra techniques.

It can be shown using K. Steiglitz and L. McBride. “A technique for theidentification of linear systems,” IEEE Transactions on AutomaticControls, vol. 10, pp. 461-464, 1965, which is incorporated herein byreference, that the objective function L_(i,init) ^((s)) can be viewedas a suitable approximation to the nonlinear LS objective function thatfollows from the relation in equation (94). Also, it should be notedthat k_(i,2) ⁽⁰⁾+k_(i,2) ⁽⁰⁾=k _(i) ⁽⁰⁾. In our experiments, we choseJ₁=20, β⁽⁰⁾=0.95, and the scalar “factor” to be 0.1 because typicallyk_(i,2) is about 10% of k_(i,3).

3. Initial Estimates for the Parameters: c, k₁

Given the initial estimates r⁽⁰⁾, l⁽⁰⁾, k_(i,2) ⁽⁰⁾, k_(i,3) ⁽⁰⁾, wegenerate initial estimates for c and k by solving the followingminimization problem

$\begin{matrix}{\left( {c^{(0)},k_{1}^{(0)}} \right)\overset{\Delta}{=}{\arg\mspace{14mu}{\min\limits_{{c \geq 0},{k_{1} \geq 0}}{J\left( {r^{(0)},l^{(0)},c,k_{1},k_{2}^{(0)},k_{3}^{(0)}} \right)}}}} & (103)\end{matrix}$

As discussed above, this problem can be solved using the ISRA algorithmas shown below in Algorithm 2 (see also equations (61) and (62)).

Algorithm 2 Initialization Procedure for the Parameters c and k₁ RunAlgorithm 1 to get initial estimates: r⁽⁰⁾, l⁽⁰⁾, k₂ ⁽⁰⁾, and k₃ ⁽⁰⁾ Forall i = 0, 1, . . ., I - 1 do  θ_(i) ⁽⁰⁾ ← [1,1,1]^(T)  Get A_(i) ⁽⁰⁾using equation (57)  For all j = 1, 2, . . ., J₂ - 1   $\left. \theta_{i,k}^{(0)}\leftarrow{\theta_{i,k}^{(0)}\frac{\left\lbrack {\left( A_{i}^{(0)} \right)^{T}d_{i}} \right\rbrack_{k}}{\left\lbrack {\left( A_{i}^{(0)} \right)^{T}\left( {A_{i}^{(0)}{\theta_{i}}^{(0)}} \right)} \right\rbrack_{k}}} \right.,{k = 1},2,3$ end for end forIn our experiments, we used J₂=500.

F. Experimental Results

We tested the IV-MBF algorithm, which is summarized in Algorithm 3below, using real patient data provided by Dr. John Votaw, Professor andVice Chair for Research, Director of Physics and Computing in Radiology,and Professor of Radiology and Physics at the Emory University School ofMedicine.

Algorithm 3 IV-MBF Algorithm Run Algorithm 1 and Algorithm 2 to getinitial estimates for r⁽⁰⁾, l⁽⁰⁾, c⁽⁰⁾, k₁ ⁽⁰⁾, k₂ ⁽⁰⁾, and k₃ ⁽⁰⁾ forall j = 0, 1, . . ., J₃ − 1 do  for all n = 0, 1, . . ., n_(0,r) dor^((m+1))[n] ← r⁽⁰⁾[n]  end for  for all n = 0, 1, . . ., n_(0,l) do  $\left. {l^{({m + 1})}\lbrack n\rbrack}\leftarrow{{l^{(m)}\lbrack n\rbrack}\frac{\sum\limits_{i = 1}^{I}{{d_{i}\lbrack n\rbrack}*{h\left\lbrack {{- n};\theta_{i,h}^{(m)}} \right\rbrack}}}{\sum\limits_{i = 1}^{I}{{{\overset{\_}{d}}_{i}\left\lbrack {{n;s_{r}^{(m)}},s_{l}^{(m)},c_{i,l}^{(m)},\theta_{i,h}^{(m)}} \right\rbrack}*{h\left\lbrack {{- n};\theta_{i,h}^{(m)}} \right\rbrack}}}} \right.,$ end for  Compute β^((m+1)) using (74) and one-dimensional line search for all n = n_(0,r) + 1, n_(0,r) + 2, ..., N − 1 do r^((m+1))[n] ←r⁽⁰⁾[n_(0,r)](β^((m+1)))^((n−n) ^(0,r) ⁾  end for  for all n = n_(0,l) +1, n_(0,l) + 2, ..., N − 1 do l^((m+1))[n] ←l^((m+1))[n_(0,l)](β^((m+1)))^((n−n) ^(0,l) ⁾  end for  for all I = 1,2, ..., I do Compute k_(i,2) ^((m+1)) and k_(i,3) ^((m+1)) usingequations (47) and (48), and a one-dimensional line search$\left. c_{i,k}^{({m + 1})}\leftarrow{c_{i,k}^{(m)}\frac{\left\lbrack {\left( A_{i}^{({m + 1})} \right)^{T}d_{i}} \right\rbrack_{k}}{\left\lbrack {\left( A_{i}^{({m + 1})} \right)^{T}\left( {A_{i}^{({m + 1})}\theta_{i}^{(m)}} \right)} \right\rbrack_{k}}} \right.,{k = 1},2$$\left. k_{i,1}^{({m + 1})}\leftarrow\frac{{k_{i,1}^{(m)}\left\lbrack {\left( A_{i}^{({m + 1})} \right)^{T}d_{i}} \right\rbrack}_{3} + {2\lambda\; k_{i,1}^{(m)}{\Sigma_{j \in N_{i}}\left( {k_{i,1}^{(m)} + k_{j,1}^{(m)}} \right)}}}{\left\lbrack {\left( A_{i}^{({m + 1})} \right)^{T}\left( {A_{i}^{({m + 1})}\theta_{i}^{(m)}} \right)} \right\rbrack_{3} + {4\lambda\; k_{i,1}^{(m)}{N_{i}}}} \right.$ end for end for

An alternative implementation simply uses the initial estimates r⁽⁰⁾,l⁽⁰⁾ after subtracting off activity due to the myocardium, as the finalRV and LV tissue curves, respectively.

1. Description of Experimental Results Analysis

The dynamic PET data set comes from scanning an unhealthy patient atrest using a standard protocol consisting of 20 scans of duration 3 sec,followed by 5 scans of duration 12 sec, and ending with 6 scans ofduration 30 sec (i.e., 31 sub-scans in total). The scanner used in thestudy has 22 planes and the reconstructed cardiac images are of size42×30. Therefore, the data set consists of 31 images of size 42×30 perplane.

Let a_(i,j) denote the measured activity in the i^(th) myocardium voxelduring the j^(th) sub-scan. Then, {a_(i,1), a_(i,2), . . . , a_(i,J)}for J=31 is the i^(th) myocardium TAC in the dynamic PET data set. Themeasured {a_(i,j)} can be modeled as

$\begin{matrix}{a_{i,j} = \left\{ \begin{matrix}{{\frac{1}{T_{1}}{\int_{u_{j,1} - T_{1}}^{u_{j,1}}{{a_{i}(t)}\ d\; t}}},{j = 1},2,\ldots,20} \\{{\frac{1}{T_{2}}{\int_{u_{j,2} - T_{2}}^{u_{j,2}}{{a_{i}(t)}\ d\; t}}},{j = 21},22,\ldots,25} \\{{\frac{1}{T_{3}}{\int_{u_{j,3} - T_{3}}^{u_{j,3}}{{a_{i}(t)}\ d\; t}}},{j = 26},27,\ldots,31}\end{matrix} \right.} & (103)\end{matrix}$

where a_(i)(t) is the i^(th) ideal continuous-time myocardium TAC (see(6)), and T₁=3 sec, T₂=12 sec, and T₃=30 sec are the sub-scan durations,and u_(j,1)

jT₁, u_(j,2)

20T₁+(j−20) T₂, and u_(j,3)

20T₁+5T₂+(j−25) T₃. The IV-MBF algorithm is based on regularly sampledmyocardium TAC data (i.e., d_(i)[n]=a_(i)(nT)). Consequently, themeasured myocardium data, {a_(i,j)}, must be pre-processed in order toestimate the regularly sampled myocardium TAC data {d_(i)[n]}. A popularapproach is to first assume the measured myocardium TAC data is nearlypiece-wise linear. Under this assumption and as discussed in W. G.Kuhle, G. Porenta, S. C. Huang, D. Buxton. S. S. Gambhir, H. Hansen, M.E. Phelps, and H. R. Schelbert, “Quantification of regional myocardialblood flow using 13N-ammonia and reoriented dynamic positron emissiontomographic imaging,” Circulation, vol. 86, pp. 1004-17, 1992, which isincorporated by reference, it follows that the values of a_(i)(t) at themidpoints of the sub-scan windows are approximately equal toa _(i)(u _(j,1) −T ₁/2)≈a _(i,j) , j=1, 2, . . . 20  (104)a _(i)(u _(j,2) −T ₂/2)≈a _(i,j) , j=21, 22, . . . 25  (105)a _(i)(u _(j,3) −T ₃/2)≈a _(i,j) , j=26, 27, . . . 31  (106)

Thus, the regularly sampled myocardium TAC data can be estimated fromthe measured myocardium TAC data using interpolation. Specifically, weuse the “known” values for a_(i)(t), which are{a_(i)(u_(j,1)−T₁/2)}_(j=1) ²⁰, {a_(i)(u_(j,2)−T₂/2)}_(j=21) ²⁵, and{(a_(i)(u_(j,3)−T₃/2)}_(j=26) ³¹, and linear interpolation to obtainestimates for the regularly sampled myocardium TACs. Note, in ourexperiment we used T_(s)=0.5 sec for the sampling interval. It should bementioned that the approach described above for generating regularlysampled myocardium TACs from the measured myocardium TAC data would alsobe used to generate the regularly sampled RV and LV TACs.

In the initialization procedure, the initial sampled RV and LV tissuecurves were obtained using equations (92) and (93), and the modifiedSteiglitz-McBride algorithm was run for J₁=20 iterations to get initialestimates for k₂ and k₃ (see Algorithm 1). Also, J₂=500 iterations ofthe ISRA algorithm were used to determine initial estimates for c and K₁(see Algorithm 2). The IV-MBF algorithm was run for J₃=300 iterationsand all one-dimensional line searches were performed using a variant ofthe golden section method. Finally, the penalty parameter was chosenexperimentally to be λ=100.

2. Discussion of FIGS.

Some results for the IV-MBF algorithm are shown in FIGS. 3 and 4. InFIG. 3, the MBF estimates for the 9^(th) plane are displayed as an imagein FIG. 3A while the corresponding myocardium mask, which specifies themyocardium voxels, is shown in FIG. 3B. Also, the PET image that resultsfrom averaging the images for the last six sub-scans and the histogramof the MBF estimates are shown in FIGS. 3C and 3D, respectively. Thesame information for plane 16 is provided in FIG. 4. Because the data isfrom a patient the true MBF values are not known. However, Dr. Votawreviewed the results and considered them to be promising and noted thatthe range of the MBF estimates was consistent with his expectation andhuman physiology.

For a specific myocardium voxel, FIG. 5A is a plot of the measuredinterpolated TAC (solid) and the corresponding estimated TAC using theIV-MBF algorithm (dashed-dotted). Also, the signal components due toactivity in the RV (dashed) and LV (dotted), and activity in the free(solid) and trapped (dashed-dotted) states are shown in FIG. 5B. Lastly,a plot of the LS objective function, which decreases monotonically asexpected, is shown in FIG. 5C.

G. Proof of Result 1

Let z₁[n], z₂[n], and w[n] be positive, casual sequences of length K.The term (z₁[n]+z₂[n]*w[n])² can be expressed as

$\begin{matrix}{\left( {{z_{1}\lbrack n\rbrack} + {{z_{2}\lbrack n\rbrack}*{w\lbrack n\rbrack}}} \right)^{2} = \left\{ {{\Delta^{old}\lbrack n\rbrack}\left( {{\frac{z_{1}^{old}\lbrack n\rbrack}{\Delta^{old}\lbrack n\rbrack}\frac{z_{1}\lbrack n\rbrack}{z_{1}^{old}\lbrack n\rbrack}} + {\frac{{z_{2}^{old}\lbrack n\rbrack}*{w\lbrack n\rbrack}}{\Delta^{old}\lbrack n\rbrack}\frac{{z_{n}\lbrack n\rbrack}*{w\lbrack n\rbrack}}{{z_{2}^{old}\lbrack n\rbrack}*{w\lbrack n\rbrack}}}} \right)} \right\}^{2}} & (107)\end{matrix}$where Δ^(old)[n]

z₁ ^(old)[n]+z₂ ^(old)[n]*w[n]), and z₁ ^(old)[n] and z₂ ^(old)[n] arenonnegative, casual sequences. Due to the convex combination on theright hand side of equation (107) and the convexity of the squarefunction, it follows that

$\begin{matrix}{\left( {{z_{1}\lbrack n\rbrack} + {{z_{2}\lbrack n\rbrack}*{w\lbrack n\rbrack}}} \right)^{2} \leq {{\frac{z_{1}^{old}\lbrack n\rbrack}{\Delta^{old}\lbrack n\rbrack}\left( \frac{{z_{1}\lbrack n\rbrack}{\Delta^{({old})}\lbrack n\rbrack}}{z_{1}^{old}\lbrack n\rbrack} \right)^{2}} + {\frac{{z_{2}^{old}\lbrack n\rbrack}*{w\lbrack n\rbrack}}{\Delta^{old}\lbrack n\rbrack}\left( \frac{\left( {{z_{2}\lbrack n\rbrack}*{w\lbrack n\rbrack}} \right){\Delta^{old}\lbrack n\rbrack}}{{z_{2}^{old}\lbrack n\rbrack}*{w\lbrack n\rbrack}} \right)^{2}}}} & (108) \\{\leq {{\Delta^{({old})}\lbrack n\rbrack}\left\{ {\frac{\left( {z_{1}\lbrack n\rbrack} \right)^{2}}{z_{1}^{old}\lbrack n\rbrack} + \frac{\left( {{z_{2}\lbrack n\rbrack}*{w\lbrack n\rbrack}} \right)^{2}}{{z_{2}^{old}\lbrack n\rbrack}*{w\lbrack n\rbrack}}} \right\}}} & (109)\end{matrix}$Now, we determine a majorizing function for the sequence (z₂[n]*w[n])².

The square of the convolution of casual, nonnegative sequences z₂[n] andw[n] can be written as

$\begin{matrix}{\left( {{z_{2}\lbrack n\rbrack}*{w\lbrack n\rbrack}} \right)^{2} = \left( {\sum\limits_{k = 0}^{K - 1}\;{{z_{2}\lbrack k\rbrack}{w\left\lbrack {n - k} \right\rbrack}}} \right)^{2}} & (110) \\{= \left\{ {\sum\limits_{k = 0}^{K - 1}\;{{\gamma_{nk}\left( {{z_{2}^{old}\lbrack n\rbrack}*{w\lbrack n\rbrack}} \right)}\left( \frac{z_{2}\lbrack k\rbrack}{z_{2}^{old}\lbrack k\rbrack} \right)}} \right\}^{2}} & (111)\end{matrix}$where z^(old)[n] is a casual, nonnegative sequence and

$\begin{matrix}{\gamma_{nk}\overset{\Delta}{=}\frac{{z^{old}\lbrack k\rbrack}{w\left\lbrack {n - k} \right\rbrack}}{{z^{old}\lbrack n\rbrack}*{w\lbrack n\rbrack}}} & (112)\end{matrix}$Because γ_(nk)≥0 for all n, k, and Σ_(k=0) ^(K) γ_(nk)=1 for all n, theconvexity property of the square function can be exploited to yield thedesired result

$\begin{matrix}{\left( {{z\lbrack n\rbrack}*{w\lbrack n\rbrack}} \right)^{2} \leq {\sum\limits_{k = 0}^{K - 1}\;{\gamma_{nk}\left\{ {{z^{old}\lbrack n\rbrack}*{w\lbrack n\rbrack}\left( \frac{z\lbrack k\rbrack}{z_{2}^{old}\lbrack k\rbrack} \right)} \right\}^{2}}}} & (113) \\{\leq {\left( {{z^{old}\lbrack n\rbrack}*{w\lbrack n\rbrack}} \right){\sum\limits_{k = 0}^{K - 1}\;{{w\left\lbrack {n - k} \right\rbrack}\frac{\left( {z\lbrack k\rbrack} \right)^{2}}{z_{2}^{old}\lbrack k\rbrack}}}}} & (114)\end{matrix}$Finally, Result 1 is obtained by replacing (z[n]*w[n])² in (109) withthe majorizing function defined by (114).

VII. Implementation of the Algorithms

The methods and techniques described in the various approaches above maybe utilized, implemented, and/or run on many different types of systems,including for example computers, game consoles, entertainment systems,etc. Referring to FIG. 6, a system 400 is illustrated that may be usedfor any such implementations. One or more components of the system 400may be used for implementing any system or device mentioned above, suchas for example even a handheld device. However, the use of the system400 or any portion thereof is not necessarily required.

By way of example, the system 400 may include, but is not required toinclude, a central processing unit (CPU) 410, a random access memory(RAM) 420, and a mass storage unit 430, such as a disk drive. The system400 may be coupled to, or integrated with, any of the other componentsdescribed herein, such as an input device 450, 460 and other inputdevice 470. The system 400 comprises an example of a processor basedsystem. The CPU 410 may be used to execute or assist in executing thesteps of the methods and techniques described herein. In one approach,the system 400 may further comprise a graphics processing unit toexecute or assist in executing the steps of the methods and techniquesdescribed herein. In some embodiments, the input device 450 may comprisea first touch sensitive panel and the input device 460 may comprise asecond touch sensitive panel. Furthermore, in another aspect, the system400 comprises another input device 460 that may comprise other userinput means such as buttons, keyboard, mouse, joystick, and the like. Inanother aspect, other input device 460 may further comprise outputmeans, such as displays, sound emitters, light emitters, and the likeconfigured to provide feedback or output to a user. In one embodimentone or more of the input device 450, input device 460 and other inputdevice 470 comprise display functionality. In one embodiment variousprogram content, images, shadows, lighting, and the like may be renderedon one or more of the input device 450, 460 and other input device 470.

The mass storage unit 430 may include or comprise any type of computerreadable storage or recording medium or media. The computer readablestorage or recording medium or media may be fixed in the mass storageunit 430, or the mass storage unit 430 may optionally include externalmemory 470, such as a digital video disk (DVD). Blu-ray disc, compactdisk (CD), USB storage device, floppy disk, or other media. By way ofexample, the mass storage unit 430 may comprise a disk drive, a harddisk drive, flash memory device, USB storage device, Blu-ray disc drive,DVD drive, CD drive, floppy disk drive, and the like. The mass storageunit 430 or external memory 470 may be used for storing program code ormacros that implement the methods and techniques described herein.

Thus, external memory 470 may optionally be used with the mass storageunit 430, which may be used for storing program code that implements themethods and techniques described herein. However, any of the storagedevices, such as the RAM 420 or mass storage unit 430, may be used forstoring such program code. For example, any of such storage devices mayserve as a tangible computer readable storage medium for storing orembodying a computer program for causing a console, system, computer, orother processor based system to execute or perform the steps of any ofthe methods, code, and/or techniques described herein. Furthermore, anyof the storage devices, such as the RAM 420 or mass storage unit 430,may be used for storing any needed database(s), gestures, lists, macros,etc.

In some embodiments, one or more of the embodiments, methods,approaches, and/or techniques described above may be implemented in acomputer program executable by a processor based system. By way ofexample, such processor based system may comprise the processor basedsystem 400, or a computer, console, graphics workstation, and the like.Such computer program may be used for executing various steps and/orfeatures of the above-described methods and/or techniques. That is, thecomputer program may be adapted to cause or configure a processor basedsystem to execute and achieve the functions described above. Forexample, such computer program may be used for implementing anyembodiment of the above-described steps or techniques for performing atask at the handheld device. As another example, such computer programmay be used for implementing any type of tool or similar utility thatuses any one or more of the above described embodiments, methods,approaches, and/or techniques. In some embodiments, the computer programmay comprise a computer simulation, or system software such as anoperating system, BIOS, macro, or other utility. In some embodiments,program code macros, modules, loops, subroutines, etc., within thecomputer program may be used for executing various steps and/or featuresof the above-described methods and/or techniques. In some embodiments,the computer program may be stored or embodied on a computer readablestorage or recording medium or media, such as any of the computerreadable storage or recording medium or media described herein.

Therefore, in some embodiments a computer program product comprising anon-transitory medium for embodying a computer program for input to acomputer and a computer program embodied in the medium for causing thecomputer to perform or execute steps comprising any one or more of thesteps involved in any one or more of the embodiments, methods,approaches, and/or techniques described herein.

While the methods and systems have been described in conjunction withspecific embodiments, it is evident that many alternatives,modifications, and variations will be apparent to those skilled in theart in light of the foregoing description.

What is claimed is:
 1. A method for estimating myocardial blood flow,the method comprising: a processing device applying a pharmacologicalkinetic model to a data set stored in a storage device, the data setderived from an imaging technique based on monitoring fluid basedtracers in a left ventricle, a right ventricle, and myocardium, whereinthe pharmacological kinetic model includes: incorporating athree-compartment kinetic model of changing concentrations of boundfluid based tracers, unbound fluid based tracers, and blood plasma fluidbased tracers into a standard factor analysis of dynamic structuresmodel without assumption that a right ventricle tissue curve and a leftventricle tissue curve obey a particular mathematical relationship bydiscretizing a continuous time model of the changing concentrations ofthe bound fluid based tracers, the unbound fluid based tracers, and theblood plasma into a discrete time model using non-negative iteratevalues for the right ventricle tissue curve and the left ventricletissue curve and non-negative factor weights, and; the processing deviceoutputting a processed data set based on the application of thepharmacological kinetic model to the data set for providing arepresentation of blood flow in the myocardium.
 2. The method of claim 1further comprising the processing device applying a least squaresobjective function to obtain estimates for parameters of thepharmacological kinetic model.
 3. The method of claim 2 furthercomprising the processing device minimizing the least squares objectivefunction.
 4. The method of claim 3 further comprising the processingdevice minimizing the least squares objective function by applying amajorize-minimize optimization technique to iteratively estimate theright ventricle tissue curve and the left ventricle tissue curve incombination with estimating the blood flow in the myocardium, whereinestimating the blood flow in the myocardium using iteratively updatedestimates of the right ventricle tissue curve and the left ventricletissue curve.
 5. The method of claim 1 further comprising the processingdevice smoothing estimates for blood flow for a given voxel by applyinga limiting factor or penalty factor on data from voxels located within agiven distance from the given voxel.
 6. The method of claim 1 furthercomprising the processing device determining initial estimates for aninitial estimated right ventricle time activity curve and an initialestimated left ventricle time activity curve and using the initialestimates for the initial estimate right ventricle time activity curveand the initial estimated left ventricle time activity curve todetermine initial parameters of the pharmacological kinetic model. 7.The method of claim 1 wherein the pharmacological kinetic model includesa semi-parametric model configured to drive time activity curves for theright ventricle imaging activity from the fluid based tracers and theleft ventricle imaging activity from the fluid based tracers to zeroover time.
 8. A non-transitory computer readable medium storinginstructions that cause a processing device, in response to executingthe instructions, to perform operations comprising: the processingdevice applying a pharmacological kinetic model to a data set stored ina storage device, the data set derived from an imaging technique basedon monitoring fluid based tracers in a left ventricle, a rightventricle, and myocardium, wherein the pharmacological kinetic modelincludes: incorporating a three-compartment kinetic model of changingconcentrations of bound fluid based tracers, unbound fluid basedtracers, and blood plasma fluid based tracers into a standard factoranalysis of dynamic structures model without assumption that a rightventricle tissue curve and a left ventricle tissue curve obey aparticular mathematical relationship by discretizing a continuous timemodel of the changing concentrations of the bound fluid based tracers,the unbound fluid based tracers, and the blood plasma into a discretetime model using non-negative iterate values for the right ventricletissue curve and the left ventricle tissue curve and non-negative factorweights, and; the processing device outputting a processed data setbased on the application of the pharmacological kinetic model to thedata set for providing a representation of blood flow in the myocardium.9. The non-transitory computer readable medium of claim 8 furthercomprising instructions to cause the processing device to performoperations comprising applying a least squares objective function toobtain estimates for parameters of the pharmacological kinetic model.10. The non-transitory computer readable medium of claim 9 furthercomprising instructions to cause the processing device to performoperations comprising minimizing the least squares objective function.11. The non-transitory computer readable medium of claim 10 furthercomprising instructions to cause the processing device to performoperations comprising minimizing the least squares objective function byapplying a majorize-minimize optimization technique to iterativelyestimate the right ventricle tissue curve and the left ventricle tissuecurve in combination with estimating the blood flow in the myocardium,wherein estimating the blood flow in the myocardium using iterativelyupdated estimates of the right ventricle tissue curve and the leftventricle tissue curve.
 12. The non-transitory computer readable mediumof claim 8 further comprising instructions to cause the processingdevice to perform operations comprising smoothing estimates for bloodflow for a given voxel by applying a limiting factor or penalty factoron data from voxels located within a given distance from the givenvoxel.
 13. The non-transitory computer readable medium of claim 8further comprising instructions to cause the processing device toperform operations comprising determining initial estimates for aninitial estimated right ventricle time activity curve and an initialestimated left ventricle time activity curve and using the initialestimates for the initial estimate right ventricle time activity curveand the initial estimated left ventricle time activity curve todetermine initial parameters of the pharmacological kinetic model. 14.The non-transitory computer readable medium of claim 8 wherein thepharmacological kinetic model includes a semi-parametric modelconfigured to drive time activity curves for the right ventricle imagingactivity from the fluid based tracers and the left ventricle imagingactivity from the fluid based tracers to zero over time.
 15. Anapparatus comprising: a processing device configured to apply apharmacological kinetic model to a data set stored in a storage device,the data set derived from an imaging technique based on monitoring fluidbased tracers in a left ventricle, a right ventricle, and myocardium,wherein the pharmacological kinetic model includes: incorporating athree-compartment kinetic model of changing concentrations of boundfluid based tracers, unbound fluid based tracers, and blood plasma fluidbased tracers into a standard factor analysis of dynamic structuresmodel without assumption that a right ventricle tissue curve and a leftventricle tissue curve obey a particular mathematical relationship bydiscretizing a continuous time model of the changing concentrations ofthe bound fluid based tracers, the unbound fluid based tracers, and theblood plasma into a discrete time model using non-negative iteratevalues for the right ventricle tissue curve and the left ventricletissue curve and non-negative factor weights, and; the processing deviceconfigured to output a processed data set based on the application ofthe pharmacological kinetic model to the data set for providing arepresentation of blood flow in the myocardium.
 16. The apparatus ofclaim 15 wherein the processing device is configured to apply a leastsquares objective function to obtain estimates for parameters of thepharmacological kinetic model.
 17. The apparatus of claim 16 wherein theprocessing device is configured to minimize the least squares objectivefunction by applying a majorize-minimize optimization technique toiteratively estimate the right ventricle tissue curve and the leftventricle tissue curve in combination with estimating the blood flow inthe myocardium, wherein estimating the blood flow in the myocardiumusing iteratively updated estimates of the right ventricle tissue curveand the left ventricle tissue curve.
 18. The apparatus of claim 15wherein the processing device is configured to smooth estimates forblood flow for a given voxel by applying a limiting factor or penaltyfactor on data from voxels located within a given distance from thegiven voxel.
 19. The apparatus of claim 15 wherein the processing deviceis further configured to determine initial estimates for an initialestimated right ventricle time activity curve and an initial estimatedleft ventricle time activity curve and to use the initial estimates forthe initial estimate right ventricle time activity curve and the initialestimated left ventricle time activity curve to determine initialparameters of the pharmacological kinetic model.
 20. The apparatus ofclaim 15 wherein the pharmacological kinetic model includes asemi-parametric model configured to drive time activity curves for theright ventricle imaging activity from the fluid based tracers and theleft ventricle imaging activity from the fluid based tracers to zeroover time.